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The  research  described  in  this  final  report  was  performed  at  The  Boeing  Company, 
Saint  Louis,  Missouri,  under  Contract  F49620-96-C-001 1,  entitled  “Nonlinear  Control  of  Fighter 
Aircraft.”  The  program  was  managed  by  Dr.  Marc  Q.  Jacobs  of  the  Dynamics  and  Control 
Branch,  Directorate  of  Mathematical  and  Computer  Sciences,  Air  Force  Office  of  Scientific 
Research,  Bolling  Air  Force  Base,  DC. 

Boeing’s  program  manager  and  principal  investigator  was  Dr.  Kevin  A.  Wise.  The 
research  described  herein  was  performed  by  Dr.  Kevin  A.  Wise,  Dr.  Jackson  L.  Sedwick,  Dr. 
Yutaka  Ikeda,  Dr.  Rowena  L.  Eberhardt,  and  Mr.  Joseph  S.  Blinker. 

The  research  reported  here  was  conducted  during  the  period  June  1996  through  June 

1999. 
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Chapter  1 
Introduction 

New  aircraft  systems  are  being  designed  to  satisfy  requirements  for  a  low  radar  signature, 
low  infrared  signature,  and  low  visible  signature.  “Low  observability”  (LO)  will  be  central  to  all 
future  fighter  aircraft.  The  design  of  LO  air  combat  systems  is  a  multidisciplinary  problem  in 
aerodynamics,  control,  electromagnetics,  and  structural  design.  Critical  flight  control  research 
problems  for  these  kind  of  aircraft  exist  in  nonlinear  and  adaptive  control,  reconfigurable  control, 
multivariable  control,  performance  optimizing  control,  tailless  aircraft  control,  and  thrust 
vectoring  for  envelope  expansion.  New  control  system  design,  analysis,  and  optimization  tools 
are  needed  to  address  the  challenges  of  controlling  these  highly  nonlinear  aircraft. 

The  Nonlinear  Control  of  Fighter  Aircraft  research  was  focused  on  exploiting  recent 
developments  in  the  use  of  Linear  Matrix  Inequalities  (LMIs)  for  tailless  fighter  control  system 
design,  analysis,  and  optimization.  Specific  progress  was  made  in  the  following  three  areas: 

•  Tailless  Fighter  Flight  Control  Design 

•  Stability  Analysis  Of  Gain  Scheduled  and  Reconfigurable  Flight  Control  Systems 

•  Aeroservoelastic  Filter  Optimization 


1.1  Research  Objectives.  Accomplishments  and  Transitions 

Our  research  objectives  in  flight  control  design  addressed  the  challenges  associated  with 
tailless  fighter  system  modeling,  design,  and  analysis.  Tailless  fighters  have  aerodynamics  that 
are  very  nonlinear,  and  they  are  typically  unstable  in  multiple  axes.  Our  flight  control  design 
research  objectives  in  tailless  fighter  flight  control  system  design  address  this  nonlinear  flight 
control  design  problem  using  a  linear  parameter  varying  (LPV)  approach,  and  using  LMIs, 
design  H2  and  Hm  optimal  control  systems  within  this  context,  using  the  LPV  models  to 
represent  the  nonlinear  control  system  (including  gain  scheduling). 

The  research  objectives  in  stability  analysis  focused  on  using  LMIs  for  analyzing 
conventional  gain  scheduled  control  systems  (currently  used  in  fighter  aircraft)  as  well  as 
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reconfigurable  flight  control  systems.  In  industry,  fighter  aircraft  flight  control  systems  are 
analyzed  using  conventional  linear  analysis  methods  on  aircraft  models  representing  a  trimmed 
flight  condition.  These  analyses  do  not  address  the  gain  scheduling  and  time  varying  parameters 
(some  slowly  varying,  some  fast)  that  are  actually  present  in  the  flight  control  system 
implementation  and  aircraft  dynamics. 

In  recent  years,  reconfigurable  flight  control  systems  have  emerged  and  are  being  flight 
tested  in  research  aircraft.  In  transitioning  reconfigurable  flight  control  systems  to  production 
aircraft,  engineers  will  need  analysis  tools  to  assess  stability  margins.  By  modeling 
reconfigurable  flight  control  systems  as  a  gain  scheduled  LPV  system,  LMIs  can  be  used  to 
assess  system  stability.  The  development  of  this  capability  was  addressed  under  this  research 
program. 

The  research  objectives  in  aeroservoelastic  (ASE)  filter  design  focused  on  using  LMIs  to 
design  filter  coefficients  over  a  range  of  flight  conditions.  The  rate  and  acceleration 
measurements  used  in  flight  control  systems  are  corrupted  by  flexible  body  motion.  Filters  are 
designed  to  remove  these  signals  so  that  they  are  not  amplified  by  the  control  system  feedback. 
Typically,  notch  filters  combined  with  low  pass  filters  are  used.  These  filters  must  be  robust  to 
the  aircraft’s  varying  mass  properties,  stores  (weapons),  and  flight  conditions.  Robustness, 
typically  built  into  the  design  by  using  wide  notch  filters,  usually  results  in  large  amounts  of 
phase  lag  near  the  loop  gain  crossover  frequency.  For  tailless  fighters  that  are  unstable  in 
multiple  axes,  this  phase  lag  may  be  unacceptable.  Our  research  objectives  in  this  area  were 
focused  on  developing  tools  for  ASE  filter  design  to  improve  upon  this  time  consuming  and 
difficult  problem. 


The  following  paragraphs  briefly  summarize  the  research  objective,  accomplishments  made, 
and  transitions  on  the  technology  in  each  area.  The  LMI  software  used  to  perform  this  research 
was  transitioned  from  Stanford  University  (Stephen  Boyd)  to  Boeing.  Journal  papers  [1-3]  and 
conference  papers  [4-8]  disseminating  this  research  are  listed  at  the  end  of  this  chapter  in  the 


references. 
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Tailless  Fighter  Flight  Control  Design 

The  research  objective  was  to  develop  LMI  based  design  tools  for  tailless  fighter  flight 
control  design.  Accomplishments  include:  1)  Deriving  a  new  LQR  based  approach  for  pitch 
plane  flight  control  system  design,  integrating  the  approach  into  an  automated  design  tool,  and 
transitioning  this  tool  to  the  National  Air  Intelligence  Center  (NAIC)  in  Dayton;  2)  Developing 
design  models  and  LMI  software  for  H2  and  Hm  optimal  control  system  design.  Seventeen 
different  H2  and  //„  optimal  control  system  designs  formulations  were  implemented.  Chapter 
2  summarizes  the  LMI  design  problems.  These  tools  have  been  transitioned  to  Boeing’s 
Guidance  and  Control  Technology  ERAD  for  further  maturation  and  application  to  aircraft  flight 
control  system  design  problems. 


Stability  Analysis  Of  Gain  Scheduled  and  Reconfigurable  Flight  Control  Systems 

The  research  objective  was  to  develop  a  LMI  based  proof  of  stability  tool  for  reconfigurable 
flight  control  systems.  Accomplishments  include  development  of  LPV  based  models  for 
modeling  gain  scheduled  and  reconfigurable  flight  control  systems,  formulation  of  the  LMI  for 
stability  analysis,  and  application  of  the  approach  to  Boeing’s  Tailless  Advanced  Fighter  Aircraft 
(TAFA).  The  LMI  analysis  was  applied  to  a  gain  scheduled  control  system  and  was  used  to 
assess  system  stability  in  the  presence  of  gain  scheduling  and  time  varying  model  parameters. 
The  LMI  analysis  was  also  applied  to  a  reconfigurable  flight  control  system  designed  to 
accommodate  battle  damage,  and  was  used  to  assess  system  stability  during  reconfiguration  of 
the  flight  control  system.  These  tools  were  also  transitioned  to  Boeing’s  Guidance  and  Control 
Technology  IRAD  for  further  maturation  and  application  to  aircraft  flight  control  system  design 
problems. 


Aeroservoelastic  Filter  Optimization 

The  research  objective  was  to  develop  an  LMI  based  tool  for  optimizing  ASE  filter  designs 
to  minimize  phase  lag  at  frequencies  near  the  loop  gain  crossover  frequency  subject  to 
constraints  on  providing  a  minimum  gain  attenuation  at  the  flexible  body  modes.  This  tool  must 
also  accommodate  the  aircraft’s  varying  mass  properties,  weapons  load  out,  and  flight 
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conditions.  Our  accomplishments  include  the  preliminary  development  of  an  LMI  tool  for 
designing  these  filters,  and  application  of  the  tool  to  Boeing’s  TAFA  aircraft  model.  It  was 
found  during  this  research  that  the  LMI  problems  grew  to  extremely  large  size,  saturating  the 
memory  on  our  workstations.  This  problem  caused  limited  success  in  this  area.  Chapter  4 
discusses  the  problem  set  up  and  our  results. 

1.2  Organization  of  the  Report 

Chapter  2  presents  theory  for  H2  and  Hm  optimal  flight  control  system  design,  application 
of  these  design  methods  to  the  Boeing  TAFA  aircraft,  and  linear  simulation  results  comparing 
conventional  and  LMI  designs.  A  detailed  description  of  the  Boeing  TAFA  aircraft  is  presented 
in  this  report  in  Appendix  A,  as  well  as  linear  models  at  key  flight  conditions  (also  in  the 
appendices).  Each  Chapter  contains  at  its  end  the  references  used  in  the  chapter. 

Chapter  3  details  the  development  and  application  of  LMIs  for  stability  analysis  of  gain 
scheduled  and/or  reconfigurable  flight  control  systems. 

Chapter  4  summarizes  our  progress  in  developing  a  LMI  based  tool  for  ASE  filter 
optimization. 
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Chapter  2 

Tailless  Fighter  Flight  Control  System  Design  Using  LMIs 


2.1  Introduction 

This  chapter  presents  derivations  of  H2  and  Hx  LMIs  for  controller  design  and  applies 
tnem  to  a  tailless  aircraft  at  a  single  flight  condition.  The  theoretical  results  presented  in  this 
chapter  build  the  foundation  for  using  LMIs  in  analyzing  stability  for  a  reconfigurable  (gain 
scheduled)  control  system  that  is  presented  later  in  Chapter  3. 

In  general,  H2  and  controller  design  is  easily  accomplished  using  control  system  design 
and  analysis  packages  like  MATRIXx  and  Matlab  by  solving  Riccati  equations.  The  difficult 
problem  that  these  tools  do  not  address  is  the  stability  analysis  of  gain  scheduled  and/or 
reconfigurable  control  systems  addressed  later  in  Chapter  3.  Modeling  these  systems  using  LPV 
models  and  posing  the  stability  question  in  a  LMI  framework  gives  the  engineer  a  new  tool  in  the 
analysis  of  complicated  flight  control  systems.  The  derivation  of  controller  design  formulas 
presented  in  this  chapter  aid  in  the  set  up  and  understanding  of  the  material  presented  in  Chapter 
3.  The  controller  design  formulas  are  presented  for  LTI  systems,  then  followed  by  time  varying 
systems.  The  time  varying  models  represent  models  typically  used  for  gain  scheduled  flight 
control  systems. 

2.2  H2  Controller  Design  Using  LMIs 

This  section  contains  the  derivation  and  LMI  problems  for  minimizing  the  H2  norm 
between  an  exogenous  variable  w  and  the  regulated  variable  z  applied  to  a  linear  time  invariant 
(LTI)  system.  Included  is  the  LMI  problem  for  solving  the  standard  regulator  problem,  and  LMI 
problems  for  linear  time  varying  systems. 

2.2.1  H2  -Norm  Minimization 

Consider  the  plant  described  by 
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x  —  Ax  +  B}W+ B2u  (2  1) 

z  =  Qx+Dnw+Dnu 

where  £{w(t  )wr  (£  )}=  78  (t  -t, ) ,  and  with  the  state  feedback  control  u  =  -Kx .  The 

controller  design  problem  is  to  select  the  feedback  gain  matrix  K  to  minimize  the  H2  -norm 
between  w  and  z .  This  problem  is  described  by 

mm£jjzr  (t)z(T)cftj.  (2*2) 

Define 

P(:)=E{x(l)xT  (<)}>0 
then,  using  Eq.  (2.1) 

E{zzT}=  (Ct  - DnK)P (Ct  -DaK)T  +D,A\ ■ 

Using  zT z  =  Tr  (zzr  )  =  £  z]  =  ||z||’  we  have 

/ 

£,{zrz}=  Tr  (zzT  ) 

=  Tr(C,PC;-DuKPC; - C,PKtD, £ -DaKPKTDTu-DnD^  ) 
-.Tr(ClC,P)-Tr(ClDnKP)-Tr(PKrD^Cl)+Tr(D,1KPKTDj,)+Tr(DnD;,) 

Defining 
Y  =  KP 

and  the  slack  variables 

X  =  DnKPKTDTu  >  0 
=  DnKPP~' PKt D{2  >  0 . 

=  DnYP~'YTDTn  >0 


Then,  using  the  Schur  complement,  this  matrix  inequality  can  be  written  as 


7  DnY' 

1 - 

o 

tO 

/  O' 

fo  Djr'  1 

- 1 

o 

0  P  _ 

1 

7 

ft. 

ytd]2  p_ 

- 1 

o 

_ 1 

L  YTDTn  P\ 

>0 


DnYP~xYTD*2  DUY 

rTr\T 


YTD] S 


TD’2  p 
,.y1 

>0 


>0 


X  D12Y 


Now, 


£{II4  }=Tr(zzT) 

=  Tr  (C[qp)-2Tr  (C(DUKP)+Tr  (X)+Tr  (DuDtu  ) 
The  state  transition  <E >(t)  satisfies 
0  =  (A-B2K)0,  0(0)  =  / 
with 

t 

x(t)  =  0(t)xo+0 (r )  jV1  (x  )%W (x  )dl  . 

0 

Assuming  x0  =  0 , 

P(/)=£{x(/>r(/)} 

=  <D  (t) J  j O'1  (t  )wt  (C  )B*0~t  (£  )duKs0T  (t) 

0  0 

=  $  (t)  j  <J>_1  (x  )B1B,r< t>~r  (x  )di0T  (t) 

0 

Differentiating  yields 

P  =  (A-B2K)P  +  P(A-B2K)T  +Brf 
or 
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-P+AP-B2Y+PAT  -YtB[  +BX=  0 


For  stability,  we  have 
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P  -  AP -  PAT  +  B2Y  +  YtBt2  -BxBTlt  0 . 

The  controller  design  problem  in  Eq.  (2.2)  can  now  be  reformulated  as  the  following  LMI 
problem. 

Hi  Controller  LMI  Problem 

min  )  Tr  (C[CtP)-2Tr  (Cf  DnKP)+ Tr  (X)+  Tr  (DnD?x  )dt  (2.3) 

0 

subject  to 

P-AP-PA7  +  B2Y + YtBt2  -  BX  ^  0 


with  K  =  YP~] . 


2.2.2  LMIs  For  Standard  LQR  Problems 

A  special  case  of  the  above  H2  design  problem  is  the  familiar  LQR  problem  that  minimizes 


with  QT  =Q>0  and  RT  =R>0  for  the  LTI  plant  x  =  Ax + Bu .  For  this  special  case 
_  0  Tx' 

o  r/1  l“- 


Then, 
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C,= 


Q 


■K' 


*  Ai  - 


>  A2 — 


0 

ft/2 


ClQ=Q,  Cl Dn  =  0,  Ar,A.=o,  A^: 


0 


0  0 

0  R7lKPK‘R 


X  vdvt  pA 


X  DJ 
YTDTn  P 


"0  0 

0 

= 

0  r^2kpktrT/2 

r/2y 

0  ytrT/2 

p 

This  last  expression  implies  that 

r^2kpktrT/2  r/2y  >q 
ytrtA  p  J 

Also,  that  Tr  (X)  =  Tr  {r/2 KPKtRT/2  j.  The  following  LMI  problem  solves  for  the  feedback 
gain  matrix  K. 


LQR  LMI  Problem 

mmTr(QP)+Tr(X) 


subject  to 


-AP - PAt  +  BY+ YtBt  -W>  0 


X 

YtrtA 


r/2y 

p 


(2.4) 


with  P  =  PT ,  X  =  XT,  K  =  YP~\  and  W>  0  given. 
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2.2.3  H2  -Norm  Minimization  For  Linear  Time  Varying  Systems 

Another  variation  of  the  above  tf2-norm  minimization  problem  is  to  consider  a  linear  time 
varying  (LTV)  plant  where  the  time  variations  are  introduced  by  gain  scheduling.  Consider  the 
LTV  plant  described  by 

x  =  A(t)x+B(t)u+w  (25) 

z-Q^^x+R^^u 

with  the  state  feedback  control  u  =  -K  (t)x .  The  controller  design  problem  is  to  select  the  time 
varying  feedback  gain  matrix  K(t)  to  minimize  the  H2-norm  between  w  and  z  over  a  time 
interval  [tmin ,/max  ] .  The  following  LMI  problem  solves  for  the  feedback  gain  matrix  K(t). 

LTV  LMI  Problem 

H,^JTr(eVp)+Tr(xW'  (2'6) 

subject  to 

x(i)  Ry2(t)Y(t)]o 

YT{t)RT/l{t )  P{t) 

with  P{t)  =  PT(t),  X{t)=XT(t),  ^(r)=y(/)/>',(0>  'min  <t<tmm,and  W(t)> 0  given. 


In  addition  to  solving  this  LMI  problem,  the  control  designer  can  minimize  the  ||z||2  and  keep 
this  norm  above  a  prescribed  level.  Consider  the  LTV  plant  described  by 

x  =  A(t)x  +  B(t)u  +  w 
z  =  Qy(t)x+Ry(t)u 
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with  the  state  feedback  control  u  =  -K  (f)x .  The  controller  design  problem  is  to  select  the  time 
varying  feedback  gain  matrix  K(t)  to  minimize  the  H2  -norm  between  w  and  z  over  a  time 
interval  keeping  the  |z|2  above  f(t).  The  following  LMI  problem  solves  for  the  feedback  gam 
matrix  K(t). 

LTV  LMI  Problem  For  Prescribed  Level  Of  Jj^j|2 

min  V  (2-7) 

ywoimmi™ /max 

subject  to 

Y  (0-7(0  >0 

7™,-y  (')+r(')>0 

P(t)-A(t)P(t)-P(t)Ar(')+B(t)Y(t)+Yr(t)BT(t)-W(t)>0 

X(l)  **(f)r(ol>0 
rT(t)RT/2(t)  P(t ) 

with  P(t)=PT  (t),  X(t)=XT(t),  K(l)=Y(l)P~'(t), 
given. 


2.2.4  H2  Controller  Design  and  Simulation  Results 

In  this  section  the  H2  LMI  controller  design  problem  is  applied  to  the  Boeing  Tailless 
Advanced  Fighter  Aircraft  (TAFA)  model  described  in  Appendix  A.  A  pitch  rate  command 
flight  control  system  was  designed  using  the  LMI  standard  LQR  formulation,  and  for 
comparison,  a  robust  servomechanism  flight  control  system  [7]  was  also  designed. 


LQR  Plant  matrices 

"-4.281534e+00 

1.029299e+00 
A  = 

0.0 

0.0 


1.022817e+01 

-1.889328e+00 

0.0 

0.0 


-7.686167e+01  0.0 

5.527753e-03  2.089781e-01 

0.0  1.0 

-3.457440e+03  -9.643200e+01 


-12- 


Nonlinear  Control 
Of  Fighter  Aircraft 


0.0 

j.  00 

o.o 

_3.457440e+03_ 

The  LQR  LMI  problem  is 
mnTr(QP)+Tr(X) 

subject  to 

-AP - PAt  +  BY+YTBr  -W> 0 

X  b/2y  a 
>0 

ytrtA  p 

with  P  =  PT ,  X  =  XT ,  K  =  YP~X ,  and  W  >  0  given.  The  resulting  LMI  controller  gains  are 

Klm,  =  V1- 1 768 1 4e'0 1  -4-81 9290e-02  -9.456374e-02] . 

The  Robust  Servo  feedback  gains  are 

Krs  ~  [-2.685933e-01  -5.909616e-02  -1.102195e-0l] 

Figure  2.1  illustrates  a  step  response  comparing  the  two  designs.  Both  flight  control  designs 
produce  similar  responses. 
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2.3  H„  Controller  Design  Using  LMIs 

In  this  section  LMI  problems  are  formulated  for  minimizing  the  Hm  norm  between  an 
exogenous  variable  w  and  the  regulated  variable  z  applied  to  a  linear  time  invariant  (LTI) 
systems  and  for  linear  thne  varying  systems. 


Consider  the  plant  shown  in  Figure  2.2,  described  by 

x  =  Ax+B^w+BjU 
z  =  C]x+Dnw+D12u 
y  =  C2x+D2lw 


(2.8) 


with  controller  given  by 

xk=Akxk+Bky 
u  =  Ckxk  +  Dky  ‘ 
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Figure  2.2  Controller  Design  Model 


Substituting  for  y  in  Eq.  (2.8)  results  in 

xk  =  Akxk  +  BkC2x+BkD2lw 
u  =  Ckxk  +  DkC2x  +  DkD2]w 

Closing  the  loop  with  the  plant  model  (Eq.  (2.8))  yields 
x  =  Ax + Bx  w  +  B2  (Ckxk  +  DkC2x + DkD2iw) 

xk  =  Akxk  +  BkC2x  +  BkD2Xw  (2-9) 

z  =  Qx+ Duw+Dn  (Ckxk  +  DkC2x+ DkD2lw) 

Expressing  Eq.  (2.9)  in  matrix  form  defines  the  closed  loop  system  matrices  (Acl,Bd,Cd,Dd)  as 


To  formulate  the  controller  design  problem  within  the  LMI  framework  we  would  like  to  write 
the  closed  loop  system  matrices  such  that  the  controller  matrix  parameters  ( Ak,Bk,Ck,Dk ) 


appear  affinely  in  the  description  of  closed  loop  system.  Let 
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0= 

.A 

D, 

and  let 

’A 

O' 

Aq  — 

0 

0 

‘0 

B2 

B  = 

I 

0 

Bn  = 


C  = 


Bx 

0 


c0  =  [c,  o] 


0  I 
C2  0 


Vx2  =[0  Dx2]  A>  = 


0 

lAi 


Then 

Ad  —  Aq  +  BQC 

bcI  =  b0+bqv2X 

Ccl=C0  +  PuQC 

A/  =  Al 

The  fact  that  the  controller  parameter  matrices  appear  affinely  in  the  closed  loop  system  is  key  to 
using  an  LMI  to  solve  for  the  controller  that  yields  internal  stability  and  provides  an  4  gain 

from  w  to  z 


jzTzdx  -~rj  wTyvdx 

To  proceed,  the  Bounded  Real  Lemma  can  be  used  to  turn  the  H„  controller  design  problem  into 
a  LMI. 

Lemma  1 

Consider  a  continuous  time  transfer  function  T  (s)  (not  necessarily  minimal)  with  realization 
T(s)  =  D+C(sl-A)~'  B .  The  following  statements  are  equivalent: 

i)  ||z) + C (r/ - ^4 )"'  < Y  and  Re(4(^))<0. 
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ii)  There  exists  a  symmetric  positive  definite  solution  X  to  the  LMI 

'atx+xa  XB  Ct~ 

BtX  -yl  DT  <0. 

C  D  -yl 

For  a  proof  of  this  lemma  see  references  [1]  and  [2].  The  above  LMI  can  be  shown  to  be 
equivalent  to  the  more  familiar  Algebraic  Riccati  Inequality  by  using  the  following  lemma. 

Lemma  2 

Thematrix  ®  S  <0  ifandonlyif  R<0,  and  Q-SR~lST  <0. 

ST  R_ 

Proof  of  Lemma  2 

r  q  si  [/  si?'iire-s/j"lsr  oir  i  o 

ST  R_~  0  I  0  7. 

Now,  consider  the  LMI  stated  in  condition  ii).  Partition  it  as  follows 

'ATX  +  XA  XB  .  CT' 

BtX  -yl  '  DT 

•  ••  .  <  0 

C  D  :  -7/ 

which  is  equivalent  to 

'SX+XA  D]<Q 

Br X  -y/J  \_DT\-1 

Combining  results  in 

'ArX+X4+y'lCrC  XB+y~'CrD  1 
BTX+y~'DTC  y~'(DTD-y2l)  < 

which  is  equivalent  to 
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o(D)<y 

ATX+XA+r'CTC+y(XB+r]CTD)(y2I-DTD)(BTX+r'DTC)<0 
which  is  the  Algebraic  Riccati  Inequality  associated  with  the  Bounded  Real  Lemma. 


Now,  consider  the  description  of  the  closed  loop  system  which  is  affine  in  the  controller 
parameter  0 .  Our  approach  is  to  solve  a  LMI  problem  to  obtain  X  and  © .  To  proceed, 
consider  the  following  lemma 

Lemma  3 

Given  ¥  =  'Fr€  Rm*m ,  and  two  matrices  P  and  Q  of  column  dimension  m ,  find  a  matrix  0 
such  that 

'¥  +  PtQQ+QtQP<  0  (2-10) 

Denote  WP  and  WQ  as  any  matrices  whose  columns  form  a  basis  for  the  null  space  of  P  and  Q , 
respectively.  Then  the  inequality  in  Eq.  (2.10)  is  solvable  for  0  if  and  only  if 

Wl'PWP<  0 
W£'¥Wq<  O' 

For  a  proof  of  this  lemma  see  references  [2]  and  [3].  By  combining  Lemma  1  and  Lemma  3,  we 
obtain  the  following  theorem  which  a  states  the  necessary  and  sufficient  conditions  for  the 
existence  of  an  Hm  (suboptimal)  controller. 

Theorem  1 

Consider  the  proper,  minimal  plant  P  (s ) .  Define 

7>=[B  0  Vl] 

Q=[C  V1X  0] 
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and  let  WP  and  WQ  be  two  matrices  whose  columns  span  the  null  spaces  of  P  and  Q , 
respectively.  Then  the  Hm  controller  can  be  found  if  and  only  if  there  exists  some  positive 
definite  matrix  Xd  such  that 

K*xWp  <o 

W^xWe<0 

where 


‘AoXj+XjAl 

B0 

x~'cT 

■A-d  '-'0 

BT0 

-Y I 

a7: 

1 - 

P 

il 

Ai 

-yi 

AqXc1  +  XdAc 

XdB0 

cl' 

Vxd  = 

BT0Xd 

-Y  / 

Ari 

Co 

A, 

-y! 

Proof:  From  Lemma  1 ,  the  controller  K(s)  =  Dk+Ck  (si  -  Ak  )"'  Bk  is  an  Hm  controller  if  and 


only  if  the  LMI 


ATdXd  +  XdAd  X dBd 
BTdXd  -XI 

Ccl  A, 


cl 

Dl 


<0 


-Y I 


holds  for  Xcl<0.  Substituting foi  {Ad,Bd,Cd,Dd)  gives 


'4xd+xclA0  XdB0  CT0' 

’cTeT&xd+xdA0 

0 

cTeTvxl 

'XdBQC 

xdeec 

BT0Xd  -yl  Dl 

+ 

v2\eBTxd 

0 

v2\QTm 

+ 

0 

0 

Cq  Ai  -yJ_ 

0 

0 

0 

VX2®C 

VX2QD2X 

XdB,  Cl' 

'cT~ 

\XdB' 

Blxct  -yl  Dl 1 

+ 

VT 

^21 

QT[BTXd  0  2%]+ 

0 

Q  A.  -I1, 

0 

y - „ - * 

vxd 

.  Aa  . 

e[c  V2X  0]<0 

a 
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(2.11) 


'i>Xd+QTQTPxd+Pxle£<0 

Note  that  PXd  -\j3r  Xd  0  Zg]  where  P -\B  0  Thus, 

X  0  °' 

P  =P  0  70. 

0  0  O 

The  left  null  space  of  P  is  related  to  the  left  nullspace  of  PXd  by 

~X-J  0  o' 

W p  =  0  1  0  wp. 

PXd 

0  0  o 

Now,  using  Lemma  3,  we  can  eliminate  ©  from  Eq.  (2.1 1)  obtaining 

'xj  o  ol  \x-J  0  O' 

W}  0  I  0  '¥Xd  0  10  WP<  0 

0  0  O  0  0  O 

v  - ✓ - - — - - ‘  J 

■  <t>.v„ 

w;*XdwP<  o 

and 

ft//1' ■ 

Theorem  1  says  that  the  set  of  7/„  controllers  is  non-empty  if  and  only  if  there  exists  a  matrix 
Xc,  satisfying  WP<&XWP<  0  and  xWa  <0.  This  characterization  involves  both  Xcl  and 
it’s  inverse  X~j .  These  conditions  can  be  reduced  to  the  more  familiar  77..  Riccati  equations  as 


follows. 


Theorem  2 
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Consider  the  proper,  minimal,  plant  -P  (i1)  of  order  n.  Let  Wn  and  W2]  denote  null  spaces  of 
(l-D+2Dn)B2  and  (/ - £>21£>2+, )C2 ,  respectively.  Then,  the  Hm  controller  (for  fixed  y  )  is 
solvable  if  and  only  if 

i)  max{a(A,),o(Z)n)}<Y 

ii)  There  exists  symmetric  positive  definite  matrices  R  and  S  such  that 


Y  / 

-bTn  Y'JUr. 


(2.12) 


AS1  \ 


(2.13) 


where 

b2  =  b2d;2  a=a-b2c ,  b1  =  b1-b2du 

Cj  =  (l — Dl2Dl2)Cj  Aj  =(j~  A2Az)Ai 

c2  =  d;2c2  a=a-b1c2  A=Q-aA 
A,=Ai(2-a*a) 

and  rank(I-RS)<k  (where  £  is  the  order  of  the  controller).  Using  X  =  A-iT1  and  T=yS“\ 
along  with  the  simplifying  assumptions 

A,  =  o  A'tAi  A]=[2  o]  A,[Ar,  Ar]=[2  «] 

the  above  Algebraic  Riccati  Inequalities  (ARIs)  can  be  written  as 
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ATX  +  XA  +  X(r1B,B[ -B2B;)X+C,Tq<0 
AY  +YA7  +r(r2cfc, -Ct2C,)Y+B,B[  <0 
x>o  r>o  p (xy)<y2 


which  are  the  more  familiar  //„  Riccati  expressions. 

The  proof  of  Theorem  2  is  based  upon  Theorem  1  and  introduces  notation  and  algebra  common 
to  much  of  the  published  literature  on  LMIs. 

Proof:  From  Theorem  1  the  set  of  H„  controllers  in  non  empty  if  and  only  if  W^XWV  <  0 
and  Wl'Yy  Wo<0  for  some  Xcl  >0  of  dimension  (n+k)x(n  +  k)  (k  is  the  order  of  the 
controller).  To  express  these  inequalities  in  terms  of  the  plant  state  space  model,  partition  Xcl 
and  X~J  as 


*c/  = 


S  N 

Nt  * 


and  X;}  = 


R  M 
Mt  * 


Substituting  these  into  Q>Xd  yields 


4>, 


ar+rat  am 

MtAt  0 


b! 


0 


A 

0 

-Y I 


rc ; 


T  r'T 


MC] 

Dl 


CXR  QM  A.  -Y I  J 


with  P  = 


0  /*  0  0 
B\  0  0  A  2 


.  The  nullspace  of  P ,  ker  (P  ) ,  has  the  form 


fV„  = 


Wx  0 
0  0 


W-, 
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where  ker  (P)= 


X 

W. , 


2J 


D[2] 


Wx 

lX 


=  [0  0].  Noting  that  the  second  row  of  WP  is 


zero,  the  condition  Wp<I>x  Wv  <  0  reduces  to 


X 

o' 

T 

'ar+rat 

B\ 

RC*~ 

X 

o' 

0 

I 

B? 

-yi 

A> 

0 

I 

A 

o 

w2 

0_ 

QR 

Ai 

-7/ 

X 

0 

Re-arranging  results  in 


ker(i?)  O' 
0  / 


T 

'ar+rat 

RC[ 

A 

CtR 

-yi 

Du 

l 

DTn 

-7/ 

ker(^)  0 
0  I 


<0 


The  next  step  consists  of  expressing  ker(i?)  in  terms  of  the  plant  matrices. 


ker(R)  = 


■  X 

o ' 

'  h  o' 

X 

o' 

.-AX 2 

i 

£5* 

A  2. 

0 

/ 

where  B2  =  B2Df2 ,  then  B2Ul2  =  B2D^2Un  =  0 .  Next,  we  see  that 


K  is]  7*  “ 

L  ^12  _ 

\X  0 
0  / 

=[#-££#  D^uu] 

rs 

&T  ° 

O' 

I 

“[tf-A \rfBl  D[2Ul2] 

CN 

feT  ° 

i _ 

=[(/-b;a*/K  a%] 

r  w 
rr  12 

0 

=  [0  0] 

Substituting  for  ker(J?)  gives 
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Wu  0 

o' 

T 

~ar+rat  rctx 

a' 

- 1 

5! 

K) 

o 

0" 

■BlWn  Un 

0 

cxr  -yi 

Ai 

i 

KJ 

CJ 

K> 

0 

<0 

0  0 

~7 

B[  Ari 

-y7 

o 

o 

1 

I 

Combining  terms  results  in 


+  RAT-yB7pT)wu 

W?2RCTXUX2 

UTnC,RWn 

-yUTnUn 

A2A1 

BTxWn 

DTnUn 

-y/ 

<0 


Further  simplification  can  be  made  by  using  UX2UX2 

'w?2(aR  +  RAt  -yB2BT2)wu 

KRC\ 

KA 

c,rwX2 

-yl 

A, 

b[wX2 

A. 

-yf 

=  1-  DX2DX 2 ,  resulting  in 


<0 


which  is  equivalent  to  Eq.  (2. 12).  A  similar  proof  is  constructed  for  Eq.  (2.13). 


The  following  sections  summarize  Hm  controller  design  formulas,  followed  by  design  results 
applying  the  LMI  controller  to  the  Boeing  TAFA  aircraft. 

23.1  State  Feedback  ARE  H„  Controller  Formulas 

The  //„  controller  formulas  presented  in  this  section  solve  an  algebraic  Riccati  Equation 
(ARE)  formulas  to  yield  a  state  feedback  control  u  =  -Kx .  Consider  the  following  LTI  plant 

x  =  Ax  +  Bxw+B2u 
z  =  C,x +Duw+ Dnu 

The  control  is  obtained  by  solving  the  ARE  for  P 

atp+pa+q-prp= 0.  (2-14) 

The  feedback  gain  matrix  is 

-24- 
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k = (RuBT2+RuBl  )p+(R„DZ ;  +RI1D[l  )c, 


where 


Ru=-D'DnZ 

R„  =  D*Z)*r  +  J!1!ZX! 


B,=B,D' 

A  =  A— B2C\ 


A=a-Aa. 

U  =  I-DUD+ 

C,=C/C, 


Ai —  ^Ai 
2_I = Ai'Ai  ~y2i 
a  =  a-b,zdtucx 


e=c,r(/-A,zA'i)c, 

R  =  B2%+B,ZB[ 


23.2  H„  Controller  LMI  Problem 

The  LMI  problem  summarized  here  solves  for  the  optimal  y  and  forms  the  state  feedback 
control  u  =  -Kx .  This  problem  is  similar  to  the  LMI  controller  design  problem  posed  by 
Gahinet  and  Apkarian  in  [2].  This  problem  set  up  is  applicable  to  extensions  made  in  applying 
LMIs  to  gain  scheduled  systems. 

Consider  the  LTI  plant  described  by 

x  =  Ax+Biw+B2u  (2.16) 

z  =  C]x  +  Duw+Duu 

The  controller  design  problem  is  to  select  the  feedback  gain  matrix  K  to  minimize  the  Hx  -norm 
between  w  and  z . 

H„  Controller  LMI  Problem 

miny  (2.17) 


subject  to 
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A  —  A  “ 

-RC\  -5, 

yA  -A.  >° 

-Ar.  Y4. . 

R>  0 

+  R]2B(  +  +^izAi)Q 

d+=a_2  A=a-4 a. 

A  =  A£+  u=j-dud+ 

a=a-b2cx  c}=ucx 

Ai  =udu 

z  =D^Dn-y  I 

Consider  the  following  graphic  illustrating  the  optimization  software’s  convergence  to  the 
optimal  y . 

Phase  1 

Feasible  Solution 


*,2=-£+A.Z 

Ru=D+D+T+RnZ~'Rl2 


-AR-RA7  +yB2BT2 
-CXR 


,  Nonlinear  Control 
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The  LMI  optimization  software  begins  by  finding  a  feasible  solution  within  the  feasible  region, 
denoted  by  X  in  Figure  2.3,  and  then  converges  to  the  optimal  y^ .  As  y  approaches  Ynun  this 
often  leads  to  ill  conditioning  in  the  H„  problem,  where  very  small  numerical  changes  in  y 
produce  very  large  changes  in  the  magnitude  of  the  feedback  gains.  H„  practitioners  [4]  then 
back  off  on  the  numerical  value  for  y ,  making  it  larger,  thus  making  the  feedback  gains 
reasonable  (for  implementation)  leading  to  a  sub-optimal  control.  The  next  LMI  problem  setup 
addresses  this  fact  within  the  LMI  framework. 


2.3.3  H„  Controller  LMI  Problem  With  Lower  Bound  Constraint 

The  LMI  problem  summarized  here  solves  for  a  sub-optimal  y  along  the  optimization  path 
(shown  in  Figure  2.3)  and  forms  the  state  feedback  control  u  =  -Kx .  Consider  the  LTI  plant 
described  by 

x  =  Ax+B^w+B2u  (2  18) 

z  =  Clx+Duw+Dl2u 

The  controller  design  problem  is  to  select  the  feedback  gain  matrix  K  to  minimize  the  Hx  -norm 
between  w  and  z . 


H x  Controller  LMI  Problem 

miny 

R,f 

subject  to 


A  A  m  A  A 

-AR-RAt  +yS2fl[ 

-RC{ 

A  ** 

-By 

-QR 

A 

-A. 

>0 

-if 

~i>n 

1 

R>  0 

y-y  >0 

K  =  (RuBt2  +RnB[)R  +  (RuDl2+RM)Cl 


(2.19) 


with 
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Rl2=-D+DuZ 
Ru=D+D+r  +RuZ']Rl2 


D+  =  DJ 

B2=B2D+ 

A  =  A-B2Cl 


Bx  —  Bi -4A. 
U  =  I-DnD+ 
CX=UCX 


DU=UDU 

z-^AVA.-y2/ 


and  f  >  0  given  as  the  desired  sub-optimal  gain. 

Figure  2.4  illustrates  the  solution  from  this  LMI  problem.  In  this  LMI  problem  the  software 
minimizes  y  along  the  optimization  path  subject  to  the  constraint  that  y-y  >0. 


Phase  1 

Feasible  Solution 


Figure  2.4  Controller  LMI  y  Optimization  With  Lower  Bound  Constraint 

Here  the  optimization  software  is  constrained  in  the  minimization  of  y  to  keep  y-y  >  0 .  The 
resulting  controller  using  this  design  has  feedback  gains  that  are  much  smaller  in  size,  and  as  a 
result  does  not  amplify  sensor  noise.  We  refer  to  this  as  gain  conditioning. 


2.3.4  Methods  Used  To  Improve  Numerical  Accuracy  And  Conditioning 

This  section  discusses  an  additional  constraint  that  was  added  to  the  LMI  problem  in  order  to 
satisfy  the  LMI  constraint  in  Eq.  (2.17)  to  within  an  £/ ,  and  a  weighting  matrix  added  to 

improve  numerical  conditioning  in  solving  for  the  Hx  controller. 
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Let 

AR-RAT+yB2% 

-RC( 

A 

-A 

iM)= 

-C{R 

Ik 

-A, 

A  rp 

-A, 

yk 

The  LMI  problem  in  Eq.  (2.17)  uses  the  constraint  that  L(y,R)>0.  Figure  2.5  illustrates 
bounding  L(y,R)>0  by  a  hypersphere  of  radius  e  ,i.e.  0  <L(y,R)<sI . 


Figure  2.5  Bounding  L(y,R )  With  si 


When  using  an  ARE  (Eq.  (2.14))  to  form  the  control  gain  matrix  one  can  put  the  solution  matrix 
p  back  into  the  ARE  to  see  how  well  the  ARE  is  satisfied  (usually  the  ARE  produces  a  matrix 
that  is  not  zero).  To  improve  how  well  the  inequality  is  satisfied,  we  investigated  constraining 
L(y,R)  by  adding  0<L(y,R)<sI  asan  additional  constraint  in  the  LMI  problem.  We  found 
this  had  no  effect  on  improving  how  well  the  inequality  was  satisfied. 

In  addition  to  investigating  the  above  constraint,  we  also  examined  various  scaling  matrices 
to  eliminate  ill  conditioning  in  the  R  matrix.  The  scaling  was  added  to  the  LMI  constraint  in  Eq. 
(2.17)  as  GL(y,R)GT  >  0 .  Eigenvalues  of  the  R  matrix  were  examined  to  scale  the  matrix  to 

improve  the  condition  number.  Like  before,  we  found  that  this  had  no  effect  on  the  numerical 
accuracy  in  satisfying  the  LMI  constraints. 
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2.3.5  An  Hm  Controller  LMI  Problem  With  Gain  Conditioning 

The  LMI  problem  summarized  here  optimizes  y ,  while  at  the  same  time  incorporates  a  lower 
bound  constraint  on  y  by  using  a  bias  Ay .  (This  LMI  problem  incorporates  gain  conditioning 
in  the  LMI  problem  setup  as  discussed  previously  in  Section  2.3.3. 

Consider  the  LTI  plant  described  by 

x  =  Ax+B{w+B2u  (2  20) 

z  =  Clx  +  Dnw+Duu 

The  controller  design  problem  is  to  select  the  feedback  gain  matrix  K  to  minimize  the  -norm 
between  w  and  z. 

Hx  Controller  LMI  Problem 

miny  (2-21) 

R,K,.y 

subject  to 

-AR-RAT+yB2BT2  -RC(  -3," 

-cxR  yi„t  -A,  >o 

-Ar  -A7,  _ 

3>0 

-AR,-R.AT  +&yB2Bl  -R.Cl  -B, 

-C,K  AY/,,  -A,  >0 

-Al  Ay/., 

y  -y  =  Ay 

K  =  (RnBT2  +RnBt)Rm+(RnDTn+RnDTu)C, 
with 
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Bl-Bl-B2Dn 
U  =  I-DnD+ 

C^C/C, 

_  Ai  =  udu 
z-'=D;}Dn-y2i 

and  y  >  0  given. 

Note  that  the  feedback  gain  matrix  K  is  calculated  using  the  matrix  R„ .  This  problem  set  up 
accomplishes  the  goal  of  making  the  feedback  gain  magnitudes  reasonable.  Unfortunately,  it 
makes  the  LMI  problem  size  very  large. 

2 .3.6  Controller  LMI  Problem  For  LTV  Systems  Caused  By  Gain  Scheduling 

In  Packard  and  Becker  [5,6]  the  problem  of  designing  gain  scheduled  H„  controllers  was 
addressed.  In  this  work  the  problem  was  made  simple  by  assuming  that  a  constant  global  y 
could  be  found  that  would  apply  over  the  entire  flight  envelope.  This  is  a  poor  assumption  for 
flight  control  systems  because  the  desired  dynamics/response  (for  an  aircraft  or  missile)  greatly 
varies  with  speed  and  altitude. 

Consider  the  problem  of  scheduling  the  control  gain  matrix  with  a  parameter  t  as 
follows: 

t 

Not  Necessarily  A  Uniform  Grid 

This  requires  gridding  the  parameter  space.  At  each  grid  point  a  Hm  controller  is  designed  using 
an  optimal  control  (ymm  )  and  solution  matrix  R  good  at  each  grid  point.  In  the  implementation 
the  control  is  no  longer  optimal  because  of  Y  and  R  variations  caused  by  the  scheduling.  From 

the  LTI  case,  Eq.  (2.17),  R  =yP~ 1  with  R  =  0 .  Now, 
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R=yp-'+y£(p-') 
at N 

=  yp-'+yp-'Pp-] 

This  leads  to 

yP~]PP~l  =yP~l  -R 

v  .  (2-22) 

=  *-R-R 
Y 

The  LMI  problem  of  calculating  the  optimal  control  for  linear  time  varying  systems  requires 
including  the  variations  due  to  the  scheduling.  This  leads  to  the  following  LMI  problem: 


H„  Controller  LMI  Problem 

miny 

R,y 


subject  to 


R-lR-AR-RAT+yB2BT2  -RC( 

Y 

-C,R  Y  h,  "A. 

-Btx  -Ar,  Y/n. 


>0 


R>  0 


K  =  (rubt2  +r12b(  )R  +  (RuD(2  +RnDlx  )C, 


with 


RU=-D+DUZ 

ru  =  d+d+t+rX2z-'r(2 


D+=DJ 

b2=b2d + 

A  =  A  —  B2Cx 


A. — UDw 

z-'=A';a,-y2/ 


(2.23) 


(2.24) 


Bx=Bx-B2Du 
U  =  I-DnD+ 
C,=f/C, 
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In  Packard  [5,6]  the  y  term  in  Eq.  (2.24)  is  assumed  to  be  zero  because  y  is  global,  and  then 

basis  functions  are  used  to  approximate  R .  This  simplifies  the  problem.  In  this  research  this 
term  is  included. 

In  aircraft  and  missile  applications,  the  following  process  is  recommended  for  developing 
flight  control  designs.  First,  develop  point  designs  over  the  flight  envelope  by  gridding  the 
scheduling  parameter  space.  Next,  fit  a  surface  to  the  resulting  optimal  y  designed  at  each 
discrete  point  as  shown  in  Figure  2.6 


Figure  2.6  Surface  Fit  To  Optimal  y 


This  surface  models  the  relative  “size”  of  die  y ’s,  and  the  variations  caused  by  the  scheduling 
parameters.  Next,  scale  this  surface  by  the  scalar  parameter  \i ,  also  shown  in  Figure  2.6.  Now, 
use  polynomials  as  basis  functions  to  represent  R  and  y  as  follows  (a  is  the  scheduling 
variable): 

R  =  Rq  +  RfiL  +  Rp}  +  •  •  • + Rjx5 
R  =  i?j  +2R2cl  +--+5RfL* 
y=y0+yla+y2a2 

Y=Yi+2y2a 

The  LMI  problem  becomes 

min  y  (2.25) 


subject  to 
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R-l  R-AR-RAT+iB$ 

i 


-QR 

-Bl 


-RC( 

-a 

A 

-A, 

-A7; 

1U 

R>  0 


(2.26) 


K(t)=(RuB?  +  R12B[)R  +  (RuD[2+Rl2D[1)Cl 

A 

Y 

Now  the  LMI  optimization  variables  are  ji  and  Ri),R1,R2,"">B5.  The  term  —  in  Eq.  (2.26)  is 
now  a  known  function,  capturing  the  variations  introduced  by  the  flight  envelope. 


2.3.7  H„  Controller  Design  and  Simulation  Results 

In  this  section  a  pitch  rate  command  flight  control  system  was  designed  using  the  H„  LMI 
controller  design  problem  stated  in  Section  2.3.2.  The  design  problem  was  applied  to  the  Boeing 
Tailless  Advanced  Fighter  Aircraft  (TAFA)  model  described  in  Appendix  A.  For  comparison, 
the  standard  ARE  approach  using  y  -iteration  was  used  to  design  a  Hm  flight  control  system 
(presented  in  Section  2.3.1).  The  design  results  showed  that  both  controllers  performed  the 
same. 

Plant  matrices 


-4.281534e+00 

1.022817e+01 

-7.686167e+01 

0.0 

0.0 

0.0 

1.029299e+00 

-1.889328e+00 

5.527753e-03 

2.089781e-01 

0.0 

0.0 

0.0 

0.0 

0.0 

1.0 

0.0 

0.0 

0.0 

0.0 

-3.457440e+03 

-9.643200e+01 

0.0 

0.0 

-1. 35469  le+00 

0.0 

0.0 

0.0 

-1.275 106e-03 

0.0 

1. 50777 le+03 

0.0 

0.0 

0.0 

0.0 

4.018092e+02 
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0.0 

0.0 

0.0 

0.0 

0.0 

B2  = 

0.0 

3.457440e+03 

0.0 

0.0 

1.3 5469  le+00 

0.0 

0.0 

-9.010294e-01 

0.0 

7.557224e+01 

0.0 

0.0 

0.0 

0.0  0.0  1.503492e+00 

0.0  0.0  0.0 

-3.457440e+02  -9.643200e+00  0.0 


0.0 

1.995138e+01 

0.0 


0.0 

9.010294e-01 

A  = 

0.0 

A  = 

0.0 

3.457440e+02_ 

0.0 

The  LMI  controller  gains  are 

Klmi  ~  [*2.283433e+01  -2.242353e+01  1.840982e+O2  8.538090e+00  5.152192e+01  -5.398792e-0l] 

The  //„  ARE  feedback  gains  are 

Kare  =  [-7.916660e-01  -7.164470e-01  4.918770e-K)0  2.491234e-01  1.639036e+00  -3.1 6920 le-05] 

Figure  2.7  illustrates  a  step  response  comparing  the  two  designs.  Both  flight  control  designs 
produce  similar  responses. 
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Figure  2.7  H„  Controller  Step  Acceleration  Command  Simulation  Response 
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Chapter  3 

Stability  Analysis  Of  Reconfigurable  and  Gain  Scheduled  Flight  Control 

Systems  Using  LMIs 


3.1  Introduction 

Stability  analysis  methods  currently  used  by  industry  for  analyzing  gain  scheduled  flight 
control  systems  do  not  address  the  time  varying  parameter  dependence  of  the  models  used  in  the 
analysis.  The  proof  of  stability  for  these  conventional  control  laws  is  provided  by  linear  analyses 
(gain  and  phase  margins)  at  discrete  points  in  the  flight  envelope  and  by  extensive  numerical 
simulation  in  non-real-time  and  real-time  simulators.  However,  proving  stability  by  numerical 
simulation  may  not  find  all  the  flight  conditions  where  the  aircraft’s  stability  and  flying  qualities 
fail  to  meet  requirements.  This  approach  can  also  be  very  expensive  in  manpower  and  schedule. 
The  development  of  a  proof-of-stability  tool  capable  of  analyzing  gain  scheduled  control  laws,  as 
well  as  reconfigurable  flight  control  laws,  will  significantly  reduce  development  costs  and  will 
provide  confidence  that  the  control  laws  will  work. 

In  addition  to  providing  normal  mode  operation,  fly-by-wire  digital  flight  control  systems  are 
being  designed  to  be  reconfigurable  and  damage  adaptive  [1].  This  capability,  once  matured, 
will  greatly  improve  flight  safety.  However,  it  complicates  the  analysis  of  the  flight  control 
laws.  In  this  paper  the  reconfigurable  control  laws  are  designed  to  be  similar  to  gain  scheduled 
control  laws.  This  is  not  true  for  all  reconfigurable  control  law  designs,  and  the  analysis  methods 
developed  here  may  or  may  not  apply  to  all  reconfigurable  flight  control  systems.  In  our  case  the 
reconfigurable  control  laws  are  scheduled  with  parameters  that,  if  known,  would  stabilize  the 
aircraft  and  provide  the  pilot  with  control  over  its  trajectory. 

Conventional  linear  stability  analysis  (point  designs  with  linear  analysis)  of  reconfigurable  or 
gain  scheduled  flight  control  laws  always  raises  the  a  question  of  validity  due  to  the 
transient/time  varying  nature  of  the  problem.  For  gain  scheduled  flight  control  laws  the  gains  are 
often  scheduled  with  angle-of-attack,  which  may  not  be  slowly  varying  in  modem  fighter 
aircraft.  For  reconfigurable  controls,  the  control  laws  often  change  very  quickly  with  the 
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identification  of  a  failure  or  damage.  If  these  systems  can  be  adequately  modeled  using  time- 
varying  parameter  dependent  models,  then  linear  matrix  inequalities  (LMIs)  can  be  used  to 
analyze  stability. 

There  has  been  much  recent  progress  made  in  using  LMIs  for  control  system  analysis  and 
design  [2-7].  For  the  past  few  years  the  CDC  and  ACC  conferences  both  have  had  special 
sessions  and  workshops  in  which  LMIs  have  played  a  central  role.  Doyle  [7]  appropriately  refers 
to  these  developments  as  a  “postmodern  control  theory.”  Virtually  all  aspects  of  control  are 
being  developed:  optimal  realizations;  controller  synthesis;  norm  scaling;  multiplier  synthesis, 
and  robustness  analysis/synthesis  for  real  parameter  variations.  Boyd  [4]  outlines  many  of  the 
control  synthesis  problems  which  can  be  posed  as  LMI  problems. 

The  following  table  lists  four  categories  we  have  used  to  describe  the  scheduling  parameters 
for  parameter  dependent  models: 


Category 

Scheduling  Parameters 

I 

Known 

Slow 

II 

Known 

Not  Slow 

III 

Inaccurate 

Slow 

IV 

Inaccurate 

Not  Slow 

Category  I  models  contain  systems  whose  parameters  are  known  and  do  not  vary  with  time.  This 
is  the  analysis  model  that  is  assumed  for  conventional  gain  scheduled  flight  control  systems  used 
in  industry.  Issues  for  this  category  include:  1)  Implementation  of  the  gain  schedule  (rates)  and 
the  finite  approximation  of  a  continuum  (interpolation  strategy);  2)  Accuracy  of  local  stability 
analysis  based  upon  point  design  and  analysis. 

Category  II  contains  systems  whose  parameters  are  known  and  vary  with  time.  Here  the  gain 
scheduling  parameter  space  is  gridded  and  at  each  grid  point  there  is  a  polytope  of  parameter 
rates.  If  parameter  dependent  Lyapunov  functions  are  used  for  analysis,  the  derivative  of  the 
Lyapunov  matrix  (?)  must  be  included.  This  creates  a  partial  differential  inequality  analysis 

problem,  requiring  approximation  methods  for  solution. 

Category  III  contains  systems  whose  parameters  are  not  precisely  known  or  measurable  and 
do  not  vary  with  time.  For  this  problem  one  could  set  up  a  polytopic  linear  matrix  inequality, 
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and  prove  stability  by  finding  a  Lyapunov  function  that,  when  using  the  control  gains  at  a  fixed 
grid  point  (the  centroid),  is  valid  within  a  convex  region  about  the  centroid  grid  point  as  well. 

Category  IV  contains  systems  whose  parameters  are  not  precisely  known  or  measurable  and 
vary  with  time.  This  problem  can  be  addressed  by  coupling  the  techniques  used  to  address 
Category  II  and  III  systems. 

In  this  research  we  are  trying  to  exploit  the  recent  developments  using  LMIs  for  parameter 
dependent  systems  to  produce  a  method/tool  that  can  be  used  to  analyze  the  stability 
characteristics  of  reconfigurable  flight  control  systems.  Our  approach  requires  modeling  the 
closed  loop  system  dynamics  in  such  a  way  that  LMIs  can  be  applied.  This  approach  is  very 
similar  to  the  ongoing  research  of  Packard,  Becker,  Balas,  Gahinet,  Apkarian  [2,3]  and  others 
who  are  developing  algorithms  for  linear  parameter  varying  gain  scheduled  control  system 
design.  To  proceed  we  must  develop  a  model  of  the  aircraft  that  is  parameterized  with  the 
amount  of  battle  damage  sustained.  A  coupled  pitch-roll-yaw  model  is  developed  in  Section  2 
that  is  parameterized  with  wing  damage.  In  Section  3,  a  gain  scheduled  control  law  for  the  pitch- 
plane  is  designed  and  scheduled  with  the  damage  parameter.  In  Section  4,  our  approach  for 
analyzing  stability  using  LMIs  is  presented.  Two  tutorial  examples  using  linear  time  invariant 
systems  are  presented  illustrating  the  stability  guarantees  from  the  LMI  analysis.  The  approach 
is  then  applied  to  the  tailless  fighter  model  using  a  control  that  is  gain  scheduled  with  the  amount 
of  wing  damage.  To  reduce  the  computational  burden  associated  with  the  LMI  problems  to  be 
solved,  only  the  pitch-plane  is  analyzed.  The  final  section  discusses  conclusions. 

3.2  Tailless  Fighter  Model  With  Battle  Damage 

High  fidelity  nonlinear  simulations  are  often  used  to  produce  linear  models  of  the  aircraft’s 
dynamics.  Once  in  linear  form,  control  laws  are  designed  and  analyzed,  and  are  then 
“scheduled”  with  parameters  that  describe  the  operating  point  or  flight  condition.  The  idea  here 
is  to  create  a  set  of  linear  models  for  the  aircraft’s  dynamics  and  control  laws,  based  upon  these 
scheduling  parameters.  The  control  laws  will  then  be  designed  in  a  conventional  way,  based 
upon  point-wise  models,  but  the  set  (using  a  closed  loop  system)  will  be  used  for  analysis. 

Consider  the  linear  differential  inclusion  (LDI)  [4]  given  by 
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x  =  Qx  x(0)  =  x0  (3.1) 

where  Q  £  R"™ .  The  LDI  can  be  thought  of  as  describing  a  family  of  linear  time  varying  (LTV) 
systems.  Specific  LDIs  of  interest  include  linear  time  invariant  (LTI)  systems,  where  Q  is  a 
singleton.  For  LTI  systems,  the  LDI  reduces  to 

x  =  Ax  +  Buii+Bww 
z  =  Cx  +  Dzvu  +  Dzww 


with 


Q  = 


l\A  B“ 

l|c  A. 


When  a  is  a  polytope,  the  LDI  is  a  polytopic  LDI  (PLDI),  with  Q  described  by  a  list  of  its 
vertices  in  the  form 


Q  =  Co 


A.J 


A* 


A 

bu,l 

BWrL 

Q 

d„,l 

This  characterization  of  the  LDI  Eq.  (3.1)  is  used  in  our  Approach.  In  Boyd  [4]  the  idea  of  using 
LDIs  to  analyze  nonlinear  systems  of  the  form 


x  =  f(x,u,w,t ) 
z  =  g(x,u,w,t) 


(3.3) 


is  discussed.  If,  for  each  (x,u,  w,t)  there  is  a  G(x,u,w,t)e  Q  such  that 
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f(x,u,w,t ) 
g(x,u,w,t)_ 


X 

=  G(x,u,w,t )  u 
w 


(3.4) 


then  certain  important  properties  about  the  nonlinear  system  can  be  inferred  from  analysis  of  the 
LDI  (such  as  the  trajectories  converging  to  zero,  etc.).  Conditions  for  the  existence  of  G  are 

/(0,0,0,0)  =  0,  g(0,0,0,0)  =  0,  and 


'£/  U_  <U_ 

\x  ?“  eQ  Vx,u,w,t.. 

dg_  dg_  dg_ 

_d  x  du  8  w  „ 

Aircraft  and  missile  models  (Eq.  (3.3))  are  typically  modeled  using  forms  similar  to  Eq.(3.4). 
Here  we  are  proposing  to  model  the  aircraft  dynamics  using  Eq.  (3.4),  under  battle  damage,  and 
form  a  linear  parameter  dependent  model  of  the  dynamics. 

The  aircraft  under  study  is  the  Boeing  Tailless  Advanced  Fighter  Aircraft  (TAFA)  model. 
Brinker  [  1  ]  contains  more  descriptions  of  the  baseline  TAFA  aircraft,  control  laws,  control  .  . 
effectors,  and  damage  models. 

A  high  fidelity  six  degree-of-freedom  (6DOF)  simulation  of  the  TAFA  was  constructed  in 
MATRIXx.  Using  the  “trim”  and  “linearization”  features  of  MATRIXx,  1  g  linear  models  were 
extracted  from  the  6DOF  with  varying  levels  of  b<Je  damage  (wing  missing).  The  linear 
parameter  dependent  model  presented  in  this  sec  non  is  a  model  of  the  coupled  pitch-roll-yaw 
dynamics  at  a  low  altitude  high  speed  flight  condition  (h=0  ft.,  Mach  =  0.9)  with  battle  damage 
(wing  missing). 

The  battle  (wing)  damage  was  parameterized  by  the  variable  p0  e  [0,l] ,  where  0  represents  a 

healthy  aircraft,  and  1  represents  a  severe  damage  mode  with  the  entire  right  wing  missing. 
Intermediate  values  of  pD ,  say  0. 1 ,  represents  the  outer  10%  of  the  wing  missing.  The  linear 
parameter  dependent  state  space  model  is  given  by: 

x  =  A(pD)x+B(pD)u. 
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Using  increments  of  0.1  for  pD ,  eleven  models  were  extracted  from  the  6DOF  for  different 
levels  of  wing  damage.  The  elements  of  the  A  and  B  matrices  were  then  fit  using  polynomials  in 
pD .  Linear,  cubic,  and  quintic  polynomials  were  used  to  model  the  matrix  elements.  For  the 
quintic  polynomial,  the  matrices  are  reconstructed  as  follows: 

A  =  Af)+AipD+A2p2D+A3p3D+A4p*D+A5piD 

B  =  B0+BlpD+B2p2D  +  B3p3D+BAp*D+B5psD 

To  determine  which  polynomial  model  to  use,  the  polynomial  approximations  were  compared 
with  the  matrices  obtained  from  the  6DOF.  Error  matrices  were  formed  for  each  polynomial 
model  using  E  =  A  -  A ,  and  the  norm  of  the  error  matrix  was  calculated.  The  ||£| ,  ||£||2 ,  |£|L , 
and  |£|  norms  were  all  used  to  measure  the  accuracy.  Figures  3.1  and  3.2  show  numerical 
results  comparing  linear,  cubic,  and  quintic  polynomials  as  a  function  of  wing  damage  parameter 
Pd- 


Error  matrix  Norm 
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A  Matrix 


Figure  3.1  Polynomial  Approximation  Of  The  A  Matrix 
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B  Matrix 


Figure  3.2  Polynomial  Approximation  Of  The  B  Matrix 

Figure  3.1  shows  that  the  linear  and  cubic  polynomials  in  pD  do  not  adequately  model  the  A 
matrix  (at  pD  close  to  zero  and  at  unity).  The  quintic  polynomial  maintains  the  norm  of  the 
error  less  than  0.5,  and  was  selected  to  model  the  dynamics.  Figure  3.2  shows  that  the  wing 
damage  parameter  pD  enters  linearly  into  the  B  matrix.  As  a  result,  the  linear  polynomial 
adequately  models  the  wing  damage  control  powers.  Numerical  values  for  the  quintic 
polynomial  coefficients  can  be  obtained  electronically  by  contacting  the  first  author  at 
kevin.a.wise@boeing.com. 
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3.3  Gain  Scheduled  Flight  Control  Design 

A  pitch  rate  command  flight  control  design  was  developed  to  accommodate  the  varying 
levels  of  wing  damage  based  on  the  Robust  Servo  Linear  Quadratic  Regulator  (RSLQR) 
formulation  [8].  The  RSLQR  is  based  on  a  linear  time  invariant  (LTI)  optimal  control  solution 
that  produces  a  constant  gain  state  feedback  controller  which  is  then  gain  scheduled  with  the 
amount  of  wing  damage  p0 .  To  produce  an  output  feedback  implementation,  projective  control 
theory  [9]  is  used  to  project  the  dominant  eigenstructure  of  the  state  feedback  RSLQR  design 
into  an  output  feedback  design. 

Figure  3.3  illustrates  the  pitch  rate  command  controller  architecture.  This  controller 
architecture  uses  pitch  rate  and  normal  acceleration  feedback.  Integral  control  is  employed  on 
the  pitch  rate  command  error  to  provide  zero  steady  state  error  to  a  step  command.  An  aircraft 
with  multiple  longitudinal  plane  control  effectors  can  be  accommodated  via  a  predefined  control 
mixing,  in  which  the  mixing  can  also  vary  with  flight  condition.  This  approach  ensures  that  the 
control  effector  blending  is  smooth  across  the  flight  envelope,  and  thus  prevents  excessive 
surface  motion  due  to  gain  schedule  variations. 

The  stick  shaping  logic  required  for  implementation  is  not  designed  in  the  RSLQR  process. 
Instead  the  stick  shaping  is  set  in  the  6DOF  controller  implementation  to  achieve  a  desired  pitch 
response  sensitivity.  In  the  6DOF  implementation  the  stick  actually  commands  normal 
acceleration  which  is  then  converted  to  an  equivalent  pitch  rate  command  prior  to  entering  the 
LQR  command  structure. 
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Figure  3.3  Longitudinal  controller  structure. 


The  feedback  gains  for  this  output  feedback  controller  were  designed  using  an  automated  tool 
that  adjusts  the  LQR  penalty  matrices.  The  gain  matrix  K  for  the  output  feedback  contains  3 
gains.  It*  location  in  the  inner  loop  control  architecture  is  shown  in  Figure  3.3.  Figure  3.3  also 
shows  a  splitter  gain  used  to  feed-forward  the  command  directly  to  the  virtual  actuator.  This 
gain  has  proven  to  be  very  important  in  shaping  the  zero  dynamics  to  get  good  flying  qualities. 

Figure  3.4  shows  the  resulting  feedback  gains  and  splitter  gain  plotted  against  the  wing 
damage  parameter  pD .  The  negative  of  the  gains  were  plotted  using  a  logarithmic  scale  to 
highlight  the  magnitude  changes.  Several  of  the  gains  have  significant  magnitude  changes 
between  design  points. 
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Figure  3.4  Pitch  Controller  Output  Feedback  Gains 


Analysis  results  for  each  of  these  point  designs  is  shown  in  Figure  3.5  along  with  the  target 
values  for  the  flying  qualities.  To  evaluate  the  flying  qualities  a  low  order  equivalent  system 
(LOES)  was  fit  to  the  closed  loop  frequency  response.  The  LOES  model  is 

q  Kjs  +  L^e* 

^SlKt  S  +^9ip(0jpS+01tp 

A  maximum  liklihood  tuning  algorithm  was  used  to  adjust  the  model  parameters.  The  tuning 

algorithm  used  the  error  in  matching  the  target  parameters  (d)sp,  t,sp,  4)>  a  minimum  stability 

margin  requirement  (6  dB  and  45  deg.  phase),  and  number  of  iterations  (10)  as  stopping  criteria. 
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The  analysis  results  show  that  the  LQR  tuning  algorithm  was  able  to  obtain  identical  flying 
qualities  as  the  healthy  aircraft  for  values  of  pD  <  0.5,  but  did  not  converge  (using  a  strict  set  of 
criteria)  for  larger  values  of  pD .  At  the  larger  values  of  pD  the  phase  margin  requirement 
terminated  the  design  process.  Further  trade  studies  could  be  performed  relaxing  requirements  to 
further  improve  the  match  with  the  target  flying  qualities. 


Figure  3.5  Point  Design  Flying  Qualities  Analysis 
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3.4  Stability  Analysis  Using  LMIs 

The  key  to  using  LMIs  for  proving  stability  is  to  model  the  system  as  a  linear  parameter 
varying  (LPV)  system  and  to  use  a  standard  Lyapunov  argument  for  proving  exponential 
stability.  First,  consider  the  linear  system  given  by 

x(/)  =  ^(/)x(r).  (3-5) 

Trajectories  for  Eq.  (3.5)  can  be  expressed  using  its  state  transition  matrix  as 

x(r)  =  <I>(r,T)x(x).  (3-6) 


For  quadratic  stability  considerations,  the  Lyapunov  function  is  given  as 

V(t)=xT  {t)P{t)x(t )  (3-7) 

where  P (t )  >  0  V/  >  0 .  Substituting  Eq.  (3.6)  into  (3.7)  yields 
V  (t)  =  xT  (x  )<I>r  (t,x  )P (t )$> (t,x  )x (x  ) 


Differentiating  yields 


V(t)  =  xT  (x)Or  {t,i)[AT  {t)P{t)+P{t)A(t)+P(t)\Q>{t,x)x{x) 
=  xT  {t)[P{t)+AT  {t)P{t)+P{t)A{t)]x{t) 


If  there  exists  positive  constantse,,  e2,  and  e3  such  that  e,/  <P(t)<e2I  and 


P(t)+AT  (r)/>(/)+P(0^(0^-E3/» 


(3.9) 


then,  substituting  Eq.  (3.9)  into  Eq.  (3.8)  yields 
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V(t)<-e3xr(t)x(t)=- 


(')*(')• 


(3.10) 


Since  exI<P(t)< e2/ , 


zlXT  (t)x(t)<V(t)<e2xT  (t)x(t). 


(3.11) 


Substituting  V(t)<e2xT  (t)x(t)  into  Eq.  (3.10)  results  in 


Integrating  this  expression  yields 


InF  (/)-lnF(x)<-— (f-x) 
e2 


In 


v(t) 

1/tOj 

V(t) 

— 7-(<exp 

F(x) 


f 


-ii(r-x) 


F(r)<F(x)exp  x) 


(3.12) 


From  Eq.  (3. 1 1),  F  (x  )  <  e2xr  (x  )x (x  )  and  F  (/) >  E,xr  (t)x(t) .  Substituting  these  expressions 
into  Eq.  (3.12)  results  in 


E,xr  (t)x(t)<£2xT  (x)x(x  )exp 


f  e  ^ 


C-) 


Nonlinear  Control 
Of  Fighter  Aircraft 


Dividing  by  xT  (t  )x  (t  )  and  using  Eq.  (3.6)  we  have 


*r(0*(0 


xT(x)QT(t,  t)<E(/,t)x(t) 
xT(x)x(x) 


.E2 

<— exp 


Co 


\ 

/ 


Which  implies  that 


which  shows  that  the  trajectories  of  Eq.  (3.5)  are  bounded  exponentially. 

Like  the  previous  linear  system  in  Eq.  (3.5),  stability  for  reconfigurable  or  gain  scheduled 
control  laws  modeled  using  a  LPV  model  described  by 

x  =  i4(p)x(/)  (3.13) 

can  be  guaranteed  by  establishing  that  a  Lyapunov  jiequality  can  be  satisfied  everywhere  within 
a  region  p  e  P  of  parameter  values,  that  is  to  say  that  there  exists  a  P  >  0  such  that 

P  +  /4T(p)/>+P^(p)^0  Vp  eP 

This  analysis  however  does  not  give  any  measure  of  the  degree  of  stability  nor  the  amount  of 
stability  sacrificed  by  requiring  stability  over  a  wide  range  of  parameter  values. 

To  this  end  we  utilize  the  standard  Lyapunov  argument  which  says  that  if  there  exists  an 
e,  >  0,e2  >  0,e3  >  0  and  a  W  such  that 
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(3.14) 


6,/sirse,/. 

and 

W  +  AtW+WA<-£3I 

then  the  state  transition  matrix  d>  (t,x  )  satisfies 


Letting  y  =  and  P  =  —  this  leads  to  the  following  inequalities 

e2  £2 

p>0,  P>I,  P+ATP  +  PA<-yI 
and  with  y  >  0 ,  this  implies 


where  X  (P)  denotes  the  minimum  eigenvalue  of  P.  This  leads  to  the  following  LMI  problem 
(letting  jx  =  — y )  which  is  the  focus  of  this  paper. 


LMI  Problem 

mm  \l  (3.15) 


subject  to 
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p(  p)>o  O') 

i-p(  p)>o  OO 

^i/-p(p)-^(p)p(p)-p(p)^(p)^o  m 

for  all  p  eP.  Then  if  p<0,  the  system  is  stable  and 

Wwf*  (3'16) 

The  parameter  (i  can  be  used  as  a  measure  of  the  degree  of  stability.  The  inequality  (//)  has  been 
introduced  to  scale  the  problem,  otherwise  the  optimization  problem  would  be  homogenous  in  \l 
and  P(p).  The  following  examples  illustrate  its  use. 

Example  1 

Consider  the  linear  system  described  by  x  =  Ax  with 

-5  1 

Hi  ■*] 

Since  this  is  a  LTI  system  tH  LMI  problem  described  in  Eq.  (3.15)  does  not  contain  a  P  term. 
Also,  since  it  is  LTI  it  is  very  easy  to  determine  stability  ( \  (J)=  -4,-6 ).  However,  this  simple 

example  will  clearly  indicate  how  (I  in  Eq.  (3.16)  can  be  used  as  a  measure  of  stability,  and  the 
bound  on  the  state  transition  matrix  can  be  determined  analytically  and  computationally  (as  a 
LMI  optimization  problem).  The  constraints  in  Eq.  (3.15)  restrict  the  P  matrix  as  follows:  Let 


P\ 

Pi 

Pi 

Pi 
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From  constraint  (i)  we  have:  /,  >  0 ,  p3  >  0 ,  and  /,/3  - p\  >  0 .  From  constraint  (»)  we 
have:  1-/,  >0,  l-/3  ^0,  and  Pi)~ Pi  ^0*  Let 


to,  to 2 

TO2  TO3 


AtP+PA 


—10/?,  +2/?2 
x 


/,-10/2+A 

2/2  —  10/3 


From  the  constraint  (/«*)  we  have:  (I-to,  >0,  Ji-TOj  >0,  and  (M--to,)(p,-TOj)-to^  >0.  To 
minimize  p  by  selecting  /, ,  p2 ,  and  /3 ,  we  can  infer  from  the  symmetry  that  /,  =/3  •  Then, 

(!-/,)(!-  /3 )-  /'  =  (1  -  A  /  -  ^2  ^  0 


For  the  problem  of  minimizing  p,  the  above  inequality  is  an  equality,  reducing  to  /22  =  (1  -  /,  f . 
Using  /,  >0  and  1-/,  >0  we  have  /2  =1-/,.  Thus  minimizing  p  also  leads  to 
(p  -  to,  )(p  -  TOj )-  to^  =  0 .  Substituting  yields 


p2  -(4-24/,)p+(-96+192/,)  =  0 

which  has  a  minimum  at  p  =  -  8  no  matter  what  /,  is.  Using  0</,  <1,  /2  =  1  -  /, ,  and 
/3= /, ,  the  eigenvalues  of  P  are  at  \  =  2/,  - 1  and  ^  =  1 .  The  minimum  eigenvalue  of  P  gives 
the  tightest  bound  ( /,  =1),  P  =  /,  which  gives 


(3.18) 


(note  that  -4  is  the  right  most  \  (A)).  Figure  3.6  (a)  plots  the  actual  norm  of  the  state  transition 
matrix  (solid)  against  the  LMI  bound  in  Eq.  (3.18).  For  this  example  they  are  identical. 
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Example  2  Bound  On  State  Transition  Matrix 


Figure  3.6  LMI  bounds  on  the  state  transition  matrix. 


Example  2 

Consider  the  linear  system  described  by  x- Ax  with 


0  1 


which  has  eigenvalues  at  \  =  -\±2j .  Using  the  same  procedure  as  in  Example  1 ,  the  P matrix 
in  Fq.  <3.17)  is  calculated  as 


To  minimize  ji  in  the  LMI  (Eq.  (3.14)),  p2  =  (2-%/2)/l0.  This  yields  a  \l  =  -0.5858, 
and  A,  (P)  si,  3-2  .  The  X(P)  =  0.17157,  which  gives 
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I  e^>  =  2 4e M  P-19) 

W) 

Figure  3.6  (b)  plots  the  norm  of  the  state  transition  matrix  (solid)  against  the  LMI  bound  in  Eq. 
(3.19).  Note  for  this  example  the  bound  is  somewhat  conservative. 

3.5  Tail'.ess  Fighter  Analysis  Results 

The  LMI  problem  described  in  Eq.  (3.15)  was  used  to  analyze  the  tailless  fighter 
reconfigurable  /  gain  scheduled  flight  control  system.  The  scheduling  parameter  used  was  the 
wing  damage  parameter  pD.  Results  are  presented  using  the  category  definitions  described  in 
the  Introduction. 

The  LMI  convex  optimization  produces  the  optimal  jxand  P  matrix  for  Eq.  (3.15)  at  the 
polytope  vertices  used  in  defining  the  analysis  problem.  The  proof  of  stability  uses  the  |1  and  the 
minimum  eigenvalue  of  P ,  X  (P) ,  to  bound  the  state  transition.  Since  the  scheduling  problem  is 

infinite  dimensional  and  gridding  of  the  parameter  space  is  used  to  make  the  problem  finite 
dimensional,  one  must  still  check  how  well  the  constraints  in  Eq.  (3.15)  are  satisfied  between  the 
grid  points. 

Stability  Analysis  For  Category  I  Systems 

Category  I  systems  are  linear  systems  whose  scheduling  parameters  are  known  and  do  not 
vary  with  time.  A  conventional  linear  analysis  on  each  discrete  flight  condition  is  valid  and  is 
typically  used  to  analyze  stability. 

The  LMI  analysis  method  in  Eq.  (3.15)  was  applied  to  the  pitch-plane  tailless  fighter  system 
described  in  Section  3  at  each  discrete  value  of  pD ,  assuming  that  pD  =  0.  Here  each  discrete 
flight  condition  is  analyzed  separately.  Solving  the  1 1  LMI  problems  produces  the  following 


results: 
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Pd 

max(Re(\(A(p jf)) 

ji  From  Eq.  (4.21) 

h(P) 

0.0 

-2.4972 

-0.04161 

2.2573  x  10-4 

0.1 

-2.4048 

-0.03984 

2.1643  x  10-4 

0.2 

-2.2151 

-0.03845 

2.0933  x  10-4 

0.3 

-1.9608 

-0.03708 

2.0222  x  10-4 

0.4 

-1.6196 

-0.03577 

1.0528  x  10-* 

0.5 

-0.9905 

-0.03466 

1.8897  x  10^ 

0.6 

-0.2471 

-0.03258 

1.7671  x  10-4 

0.7 

-0.0918 

-0.02940 

1.5978  x  10*4 

0.8 

-0.0584 

-0.02616 

1.4250  x  10"4 

0.9 

-0.0447 

-0.02294 

1.2523  x  10-4 

1.0 

-0.0393 

-0.01991 

1.0893  x  10"* 

The  first  column  is  the  wing  damage  parameter.  The  second  column  is  the  maximum  real  part  of 
the  eigenvalues  for  the  closed  loop  system  matrix  using  the  gain  scheduled  control.  The  third 
and  fourth  columns  form  the  exponential  bound  described  in  Eq.  (3.16).  The  LMI  analysis 
shows  that  the  system  is  stable,  but  conservatively  bounds  the  state  transition  matrix.  Here  LMI 
analysis  is  not  really  necessary  since  the  examination  of  the  eigenvalues  of  the  closed  loop 
system  matrix  suffices  to  prove  stability. 


Stability  Analysis  For  Category  II  Systems 

Category  II  systems  are  linear  systems  whose  scheduling  parameters  are  known  and  vary 
with  time.  A  conventional  linear  analysis  at  each  discrete  flight  design  point  ignores  the  fact  that 
the  scheduling  parameters  can  vary  with  time.  Here,  polytope  vertices  are  introduced  to 
accommodate  the  time  varying  scheduling  parameter. 

Rates  of  change  of  parameters  will  produce  a  P  term  in  the  Lyapunov  inequality  in  the  LMI 
(Eq.  (3.15))  which  can  be  described  by 


•  ■r-i  d  P  . 

peP 


If  we  assume  that  the  rates  of  the  parameter  variations  p  are  known,  the  problem  becomes 
determining  the  partial  derivatives  d  P/d p  .  The  presence  of  these  terms  creates  a  linear  partial 
differential  inequality.  In  the  LMI  problem  we  must  find  a  |i  valid  for  the  entire  range  of 
parameter  values  0  <  p  <  1 .  To  solve  this  problem  we  use  basis  functions  (similar  to  [5,6])  in 

which  we  model  P(p  )  using  polynomials.  The  polynomial  model  is 
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P(  P)=|/!P' 

1=0 

where  DP  is  the  order  of  the  polynomial.  Then, 


Polytope  vertices  are  then  formed  at  the  min  and  max  rates  for  pD.  Since  we  can  not  add  the 
wing  back  on,  the  min  rate  is  modeled  as  p0  =  0 .  The  maximum  rate  was  varied  to  determine 
the  impact  on  the  LMI  analysis. 

Figure  3.7  shows  the  LMI  stability  analysis  results  for  the  Category  II  system  model.  The 
four  curves  represent  results  using  different  max  rates  (pmax  =  0.02,  1 0)  and  different  order 

polynomials  for  P  (  p  ) ,  (Dp  =5,  3 ) .  The  optimal  p  determined  from  the  convex  optimization 

was  p  =  -0.01991,  and  was  the  same  for  all  four  cases.  The  plots  in  Figure  3.7  show  how  well 
the  optimal  P  satisfies  the  LMI  constraints  at  points  in  between  the  0.1  pc  grid  spacing.  The 
first  plot  shows  the  (P  >  0)  constraint  was  satisfied.  The  second  plot  shows  that  the  (/ -  P  >  0) 
constraint  is  violated  in  betweeu  0D  =  0.9  and  1.0.  The  third  constraint  in  the  LMI, 
p/-P(p)-jfr(p)/>(p)-7>(p)J4(p)^0,  was  examined  atthe  vertices  where  p  =0  and 

where  p  =  pmax .  The  curves  show  X(A)  where  A  =  -/>(p)-^4r(p)/>(p)-i>(p)^(p),  which 
we  would  like  to  be  positive,  and  show  that  this  is  also  violated  in  between  pD  =  0.9  and  1 .0. 
However,  as  long  as  the  optimal  fi  is  big  enough  to  make 

^/_/>(p)_^(p)P(p)_Jp(p)^(p)>0  we  have  a  good  solution  to  the  LMI  problem.  The 
difference  between  p  and  X  (A)  can  be  thought  of  as  a  “stability  margin”  in  the  analysis.  Close 
examination  of  the  most  negative  point  in  Figure  3.7  (lower  right  plot)  shows  that  the  optimal  p 
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=  -0.01991  is  not  big  enough  to  cover  X(A)  =  -0.02273 .  This  result  indicates  that  the  P  matrix 
does  not  prove  stability  at  this  value  of  pD .  (Later  we  will  show  how  to  correct  this  problem.) 


Figure  3.7  Category  II  stability  analysis  results. 
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Stability  Analysis  For  Category  in  Systems 

Category  III  systems  are  linear  systems  whose  scheduling  parameters  are  not  precisely 
known  and  are  slowly  varying  with  time.  Since  the  scheduling  parameters  are  not  accurately 
known,  the  feedback  gains  that  are  scheduled  with  these  parameters  may  not  match  the  dynamics 
as  intended.  In  this  case  we  would  like  the  gains  to  provide  stability  and  performance  in  a 
convex  region  about  the  nominal  or  centroid  grid  point 

This  convex  region  must  be  formed  based  upon  estimates  of  the  accuracy  in  estimating  the 
scheduling  parameters.  For  example,  consider  the  following  scheduling  parameter  centered  at 

Pnom 


Centroid 

I  Grid  Points 

♦ 


Uncertain  Region 

Feedback  gains  have  been  designed  at  each  grid  point.  Estimates  on  the  accuracy  of  the 
scheduling  parameter  define  the  uncertain  region,  which  in  this  example  we  have  assumed  is 
smaller  than  the  grid  size.  To  analyze  stability  using  the  gains  designed  at  pnom ,  a  convex 

polytope  covering  the  uncertainty  region  is  formed  and  analyzed. 

Consider  the  parameter  dependent  system  described  by  Eq.  (3.13).  The  LMI  problem  for  this 
system  is  infinite  dimensional  because  of  the  continuous  parameters  used  for  scheduling.  The 
problem  can  be  made  finite  dimensional  by  densely  gridding  the  p  parameter  space. 

To  include  uncertainties  in  the  parameter  p  ,  consider  a  nominal  value  with  ±3ap  variations, 
described  by 

p.  =  p -3op  <p<p+3a„  =p2 
Then, 
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p=ep,+(i-e)p2,  o<e  <1.  (3-2°) 

Consider  the  case  where  the  system  matrix  ^(p)  can  be  represented  as  ^(p)  = 

A0  +  Aip+A2p2.  Using  Eq.  (3.20),  p2  =82p2+28(l-e)p1p2+(l-e)2  p2.  But  we  can  write 

P  =e2p,+2e  (l-e)PlP^  +  (l-e)2  p2.  This  allows  us  to  describe  ^4(p)  as 

A(p)=Q2\_A<)+Alpi+A2p2'j+‘2Q  (1_®  )[^0  +  4'£l2£L+-^2p1P2]  ^ 

+  (1 -8  )2  [  4,  +  4  p2  +  ^2  P2  ] 

where  the  weights  8 2 ,  28  (1  -8  ),  and  (1  -8  )2  are  non-negative  and  sum  to  unity.  This 
represents  /l(p  )  in  barycentric  coordinates  at  the  point  p  between  p,  and  p2.  Also,  A( p)  is 
in  the  convex  hull  formed  by  the  vertices  of  the  triangle  represented  by  the  bracketed  expressions 
in  Eq.  (3.21).  This  can  be  shown  graphically  as  follows: 


This  can  be  generalized  to  any  degree  of  polynomial.  For  system  data  A(p  )  represented  by  a 

polynomial  of  degree  DA,  i.e.  A(p)='^A1pi ,  there  will  be  a  convex  polytope  that  has  DA  +1 

/=0 

vertices  that  will  enclose  the  system  data  in  the  interval  [p, ,  p2  ] .  Let  the  jth  vertex  be  denoted 
as  Vj .  Then 
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where  for  j<DA-i 


h  i 


IM 

lit! 


(pi-p2)' 


P2 


and  for  j~2.DA-  i 


*’.>Wf£TvV^ 


k+i+j-Dj 


where  the  brackets  denote  binomial  coefficients.  It  is  straight  forward  that  this  procedure  can  be 
extended  to  any  number  of  scheduling  parameters. 

Figure  3.8  shows  the  LMI  stability  analysis  results  for  category  III  systems.  Here,  the 
polytope  vertices  are  enclosed  using  the  above  method  with  3c  p  =  0.1 .  The  first  curve  listed  in 

the  legend  shows  results  from  category  II  (c  p  =  0,  Dp  =  3,  pmax  =  0.02 )  for  comparison.  The  next 

three  curves  use  (pm„  =  0.02)  and  different  order  polynomials  for  P(p),  (Dp  =5,  3) .  The 

fourth  curve  listed  incorporates  additional  constraints  to  prevent  X(A)  from  becoming  too 

negative.  By  adding  additional  LMI  constraints  in  between  0.9  and  1.0,  the  resulting  optimal  P 
matrix  better  satisfies  the  constraints.  The  cost  is  a  larger  computational  LMI  problem.  The 
optimal  \l  for  the  Category  III  model  was  |i  =  -0.01632.  We  see  from  the  figure  that  by  adding 
the  additional  constraints  in  between  pD  =  0.9  and  1.0  the  LMI  analysis  proves  stability. 
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Stability  Analysis  For  Category  IV  Systems 

Category  IV  systems  are  linear  systems  whose  scheduling  parameters  are  not  precisely  known 
and  vary  with  time.  Our  solution  approach  for  Category  IV  models  combines  the  approaches 
used  for  Category  II  and  III  systems,  forming  polytope  vertices  based  upon  the  uncertainty  in  p 
and  p  .  Figure  3.9  shows  the  LMI  stability  analysis  results.  Here,  the  polytope  vertices  were 

formed  using  3ap  =  0.1  with  (pmax  =  5),  and  different  order  polynomials  for  P(p  ), 

(DP  =1,  3) .  The  third  curve  listed  in  the  legend  added  additional  constraints  to  keep  the 
constraints  satisfied  in  between  pD  =  0.9  and  1 .0.  In  calculating  the  fourth  curve,  grid  points  at 
small  values  of  pD  were  dropped  and  new  grid  points  in  between  pD  =  0.9  and  1 .0  were  added, 
where  the  total  number  of  grid  points  was  kept  equal  to  the  original  problem  (thus  keeping  the 
size  of  the  LMI  problem  constant).  The  optimal  p  for  each  curve  listed  in  Figure  3.9  was  p  =  - 
0.01627,  -0.01586,  -0.01574,  and  -0.01632,  respectively.  The  results  prove  stability  for  the  gain 
scheduled  control  in  the  presence  of  time  varying  uncertain  parameters. 
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Figure  3.8  Category  III  stability  analysis  results. 
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3.6  Chapter  3  Conclusions 

A  parameter  dependent  model  of  a  tailless  fighter  with  battle  damage  was  developed  to 
investigate  methods  of  stability  analysis  for  reconfigurable  flight  control  systems.  Stability 
analysis  methods  using  linear  matrix  inequalities  were  presented  for  parameter  dependent  models 
characterized  by  their  scheduling  parameters.  Four  categories  of  scheduling  parameters  were 
discussed  along  with  their  respective  analysis  method.  For  reconfigurable  control  systems  that 
can  be  modeled  within  this  framework,  the  linear  matrix  inequalities  approach  used  for  stability 
analysis  provides  a  proof  of  stability  for  the  gain  scheduled  control  system. 
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Chapter  4 

Aeroservoelastic  Filter  Design  and  Optimization  Using  LMIs 


4.1  Introduction 

This  chapter  presents  the  development  of  a  LMI  based  tool  for  optimizing  aeroservoelastic 
filter  coefficients  used  in  the  filtering  inertial  measurement  signals  (angular  rates  and 
translational  accelerations).  Limited  success  was  obtained  in  developing  this  tool  due  to  the 
extremely  large  sizes  of  the  resulting  LMI  problems. 

4.2  ASE  Compensation  Filter  Design 

The  influence  of  the  structural  modes  on  the  sensor  measurements  must  be  attenuated 
through  the  use  of  aeroservoelastic  (ASE)  compensation  filters  prior  to  use  in  the  control  law, 
since  insufficient  attenuation  of  these  effects  can  lead  to  instability.  Under  the  RESTORE 
program,  linearized  models  of  the  ASE  dynamics  were  extracted  from  the  TAFA  6DOF  at  the 
high  speed  and  low  speed  flight  conditions ,  and  were  used  to  design  ASE  compensation  filters 
for  the  flight  control  system.  Under  RESTORE,  the  ASE  compensation  filters  were  designed  to 
gain  stabilize  the  system  at  frequencies  above  the  first  structural  mode  by  providing  at  least  8  dB 
of  attenuation  at  these  frequencies.  Over  design  of  these  filters  can  result  in  excessive  low 
frequency  phase  lag  which  degrades  rigid  body  stability  margins  and  flying  qualities.  As  a 
result,  the  filters  were  designed  to  meet  the  high  frequency  attenuation  goals  while  minimizing 
low  frequency  phase  lag.  The  resulting  filters  are  shown  in  Figure  4.1  and  a  corresponding 
frequency  response  in  Figure  4.2. 

Frequency  responses  of  the  linearized  ASE  dynamics  showed  the  highest  amplitude  peaks  in 
the  40-50  rad/s  region  with  slightly  lower  peaks  at  higher  frequencies.  Since  the  highest  peaks 
were  at  the  lowest  resonant  frequencies,  utilization  of  low  pass  filters  to  provide  the  required 
attenuation  was  not  a  viable  option.  A  reasonably  low  order  filter  (i.e.  <  4th  order)  would  require 
a  very  low  break  frequency,  which  would  lead  to  excessive  low  frequency  phase  lag.  A  notch 
filter  provides  a  viable  method  for  attenuating  the  low  frequency  peaks  by  concentrating  the 

attenuation  in  the  desired  frequency  range.  Keeping  the  attenuation  band  of  the  notch  filter 
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narrow  minimizes  the  associated  low  frequency  phase  lag.  Cascading  the  notch  filter  with  a  low 
pass  filter  provides  attenuation  above  the  notch  filter  frequency  to  attenuate  higher  frequency 
modes.  This  filter  combination  is  commonly  used  in  ASE  filter  design. 

The  notch  -  lag  filter  combinations  designed  for  the  TAFA  airframe  are  provided  in  Figure 
4.1.  ASE  compensation  filters  were  made  as  common  as  possible  for  all  signals  in  a  given  axis 
to  prevent  an  adverse  performance  impact  when  the  signals  are  blended.  An  example  of  this  is 
the  calculation  of  stability  axis  rates.  If  the  roll  rate  and  yaw  rate  are  processed  by  filters  with 
different  bandwidths  and  then  transformed  to  stability  axis,  the  resulting  signal  can  differ 
significantly  from  a  filtered  version  of  the  true  stability  axis  rates.  Unmatched  filters  introduce  a 
phase  discrepancy  between  the  signals  leading  to  the  degraded  results.  The  lag  filters  used  in  the 
lateral  directional  axis  were  scheduled  with  Mach  to  gain  the  required  structural  mode 
attenuation  at  low  speeds  while  preserving  good  rigid  body  stability  margins  at  high  speeds. 


04  1.0 

Mach 


Figure  4.1  TAFA  ASE  Compensation  Filters 

An  example  of  the  frequency  response  analysis  used  to  validate  the  ASE  compensation  filters 
is  provided  in  Figure  4.2.  This  data  is  for  the  longitudinal  axis  at  the  low  altitude  ingress  /  egress 
flight  condition  described  in  Reference  1. 

The  first  row  of  plots  shows  the  magnitude  of  the  frequency  response  relating  the  ASE 
contribution  to  sensed  pitch  rate  and  normal  acceleration  to  the  pitch  axis  control  effector  inputs. 
This  data  shows  high  frequency  peaks  in  the  frequency  response  curves  which  exceed  20  dB, 
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indicating  the  need  for  ASE  filtering  in  the  control  law  design.  The  highest  peaks  occur  in  the 
40-50  rad/sec  region  as  discussed  above. 

The  second  row  of  plots  shows  the  open  loop  frequency  response  with  the  loop  broken  at  the 
actuator  for  each  of  the  TAFA  pitch  axis  control  effectors.  The  open  loop  model  includes  the 
rigid  body  and  ASE  dynamics,  linearized  TAFA  baseline  control  law,  actuator  and  sensor 
dynamics,  digitization  effects,  and  time  delays.  ASE  compensation  filtering  was  not  included  in 
this  analysis  to  determine  the  need  for  filtering.  The  magnitude  frequency  responses  with  the 
loop  broken  in  the  trailing  edge  flap  (TEF)  and  pitch  thrust  vectoring  channels  show  resonant 
peaks  exceeding  0  dB  in  the  high  frequency  region  (above  20  rad/sec).  This  indicates 
insufficient  attenuation  of  the  structural  modes  and  the  potential  for  instability  (stability  depends 
on  the  phase  characteristics  in  this  frequency  band).  Since  the  phase  at  these  frequencies  is 
uncertain,  gain  stabilization  of  the  structural  modes  is  desired,  indicating  the  need  for  ASE 
compensation  filtering. 

The  last  row  of  plots  in  Figure  4.2  repeats  the  open  loop  frequency  response  analysis  with  the 
addition  of  the  ASE  compensation  filters  described  in  Figure  4.1.  These  data  show  that  these 
filters  provide  the  desired  attenuation  of  the  structural  modes,  while  maintaining  good  rigid  body 
stability  margins^  The  loop  gain  at  this  flight  condition  results  in  no  0  dB  gain  crossovers, 
thereby  yielding  an  infinite  phase  margin.  The  -180°  phase  crossovers  occur  in  the  10  -  20 
rad/sec  region,  yielding  gain  margins  in  excess  of  10  dB. 

These  ASE  compensation  filters  were  also  integrated  into  the  TAFA  6DOF  for  analysis. 
Nonlinear  simulations  with  the  ASE  dynamics  described  in  Appendix  A  was  performed  to  verify 
that  the  closed  loop  system  was  unstable  prior  to  incorporation  of  the  ASE  compensation  filters, 
and  that  the  desired  low  order  response  characteristics  (flying  qualities)  were  obtained  after  the 
filters  were  added. 


Magnitude  (dB)  Magnitude  (dB)  Pitch  Rate  (deg/s) 
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Open  Loop  Freq  Response  with  Loop  Broken  at  Actuator  Baseline  Control  Laws  without  ASE  Filtering 


Open  Loop  Freq  Response  with  Loop  Broken  at  Actuator  Baseline  Control  Laws  with  ASE  Filtering 


Figure  4.2  Open  Loop  Frequency  Response  Analysis  to  Validate  ASE  Compensation  Filter 

Design 
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4.3  Optimization  of  ASE  Filter  Coefficients 

In  the  Boeing  Phantom  Works,  a  tool  called  COMP  ASE  is  typically  used  to  optimize  ASE 
filter  coefficients  to  minimize  phase  lag  at  low  frequencies,  subject  to  satisfying  gain  attenuation 
constraints  at  the  flexible  body  modes.  The  input  to  COMP  ASE  is  a  worst-case  magnitude 
frequency  response  built  over  a  range  of  flight  conditions,  store  loadings,  and  mass  properties. 
These  filters  are  designed  to  be  robust  to  the  aircraft’s  configuration  so  that  control  law  changes 
are  not  required  when  different  weapons  are  loaded  onto  the  aircraft 

COMP  ASE  is  a  conjugate  gradient  based  optimization  tool.  It  optimizes  ASE  filter 
coefficients  to  provide  gain  attenuation  at  the  flexible  body  modes  while  minimizing  phase  lag  at 
low  frequencies.  It  is  time  consuming  to  use,  and  sometimes  has  difficulty  in  convergence. 

The  design  of  ASE  filters  are  in  the  critical  path  for  developing  an  aircraft’s  control  laws. 
Using  analytical  models  for  the  flexible  dynamics,  these  filters  are  designed  and  integrated  into 
the  control  laws.  After  completion  of  ground  vibration  testing  and  structural  mode  testing,  these 
filters  are  often  re-designed  based  on  measured  frequency  response  data.  Tool  improvements  in 
this  area  could  significantly  reduce  the  costs  and  schedule  associated  with  the  design  of  the  ASE 
filters. 

Boyd  [2,3]  investigated  using  LMIs  to  design  finite  impulse  response  filter.  Instead  of 
designing  the  frequency  response  X  (to)  of  the  filter  directly,  the  power  spectrum 

is  designed,  leading  to 

find  R  (g)  ) 

suchthat  1}  (co,  )<R((ai)<U2  (co, ),  z'  =  l,---,A/ 

J?(cq)>  0,  coe£2£[0,2jt] 

where  L  (co  )  and  U  (co  )  are  magnitude  bounds  for  the  filter,  L(co,)<  |Xo),  |  <  U  (go,  ) , 

/  =  1,  •  •  • ,  M .  This  approach  requires  M  to  be  sufficiently  large. 

The  ASE  filters  for  TAFA  (Appendix  A)  are  implemented  at  600  Hz,  which  is  the  rate  that 
the  IMU  provides  data  to  the  flight  control  system.  This  data  is  buffered  and  processed  at  100 
Hz  to  form  the  feedback  measurements  for  flight  control.  The  high  sample  rate  (600  Hz)  creates 
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a  large  number  of  frequencies,  M .  As  a  result  the  size  of  the  LMI  problem  grew  too  large  for 
numerical  optimization. 

If  the  sample  rate  for  implementation  is  lowered,  the  problem  size  does  become 
implementable.  Figure  4.3  shows  the  upper  bound  U  (co)  and  lower  bound  L(g>)  frequency 
responses  for  a  200  Hz  filter  design  using  52  frequency  points  ( M  =  52). 


Figure  4.3  200  Hz  Filter  Design  Bounds 


The  spike  at  7.32  inputs  the  notch  required  to  attenuate  the  first  bending  mode.  Using  these 
design  constraints,  the  LMI  problem  was  solved  optimizing  the  filter  tap  coefficients.  For  a  filter 
with  21  tap  coefficients  the  results  were 

x[0]=0. 789425 
x[l]=-0. 320252 
x[2]=0. 241448 
x[3]=-0. 148961 
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x[4]=0. 060187 
x[5]=-0. 027118 
x[6]=0. 026300 
x[7]=-0. 023576 
x[8]=0. 066692 
x[9]=-0. 021337 
x[10]=0. 026092 
x[ll]=0. 027752 
x[12]=0. 006417 
x[13]=0. 058977 
x[14]=0. 026731 
x[15]=0. 106196 
x[16]=0. 011379 
x[17]=0. 303340 
x[18]=-0. 077172 
x[19]=0. 555990 
x[20]=-0. 569076 


Figure  4.4  shows  the  resulting  frequency  response. 


Normalized  Frequency  (rad/s) 


Figure  4.4  200  Hz  Filter  Frequency  Response 
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Here  the  notch  filter  shown  in  Figure  4.4  would  be  combined  with  a  low  pass  filter  to  roll  off  the 
larger  gains  at  the  upper  frequencies. 

4.6  Chapter  4  Conclusions 

Our  research  in  this  area  has  shown  that  designing  FIR  filters  using  LMI  tools  is  feasible. 

Our  application  to  designing  ASE  filter  coefficients  is  probably  no  the  best  application  for  this 
theory.  The  resulting  LMI  problems  become  very  large  in  size  and  easily  exceed  the  memory  of 
most  common  work  stations.  As  a  result,  this  does  not  improve  the  toolset  for  the  control  system 
designer. 
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The  Boeing  Tailless  Advanced  Fighter  Aircraft  (TAFA),  illustrated  in  Figure  A.  1,  is  a 
conceptual  design  of  an  advanced  fighter  configuration  which  blends  an  extensive  suite  of 
conventional  and  innovative  control  effectors  to  achieve  high  agility  in  a  low  observable  design. 
The  TAFA  is  a  single  engine,  single  seat  fighter  designed  for  air  to  air  or  air  to  ground  missions. 


Figure  A.1:  Three-View  of  TAFA  Aircraft 


The  TAFA  airframe  is  characterized  by  a  chined  forebody,  symmetric  air  inlets,  and  the  lack 
of  a  vertical  tail.  The  wing  and  all  moving  canard  are  thin  and  feature  a  moderate  aft  sweep  with 
no  dihedral.  The  wing  is  equipped  with  leading  edge  passive  porosity  that  can  be  used  as  a  low 
rate  roll  control  device  during  covert  maneuvers.  The  passive  porosity  strip  consists  of  a  series 
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of  valves  which  can  be  opened  to  reduce  the  lift  on  one  wing  of  the  aircraft,  thereby  producing  a 
roll  maneuver. 

The  trailing  edge  of  the  wing  features  ailerons,  trailing  edge  flaps  and  aft  body  split  flaps. 

The  trailing  edge  flaps  provide  a  powerful  pitch  control  effector  which  can  also  be  deflected 
differentially  to  augment  the  ailerons  during  rolling  maneuvers.  If  necessary,  the  flaps  and 
ailerons  can  be  deflected  in  opposing  directions  to  generate  yawing  moments  without  inducing 
roll.  The  aft  body  split  flaps  are  “clamshell”  devices  which  consist  of  two  panels  on  each  side  of 
the  aircraft.  One  panel  opens  above  the  wing  and  the  other  below  the  wing  to  produce  yawing 
moments  while  inducing  very  little  roll.  The  all  moving  canards  are  used  as  a  low  rate  trim 
device  for  performance  optimization,  but  also  provide  supplementary  yaw  control  power  through 
differential  deflections.  In  addition,  the  canards  generate  substantial  nose  down  control 
capability  to  help  meet  control,  margin  requirements  at  high  angles  of  attack. 

The  TAFA  is  powered  by  a  moderate  bypass  ratio  turbofan  engine  equipped  with  axi- 
symmetric  thrust  vectoring.  The  pitch  and  yaw  thrust  vectoring  enhance  maneuvering 
capabilities  and  stability  augmentation.  In  addition,  compressor  bleed  air  is  routed  to  forebody 
ports  for  pneumatic  control.  The  ports  are  mounted  to  serve  as  a  yaw  control  device. 

The  TAFA  is  designed  for  fly-by-wire  control  using  hydraulically  and  electromechanically 
actuated  control  effectors.  Movable  surfaces  utilize  hydraulic  actuators  while  the  pneumatic 
forebody  blowing  and  leading  edge  passive  porosity  employ  electromechanical  actuators.  The 
system  is  designed  to  be  tri-redundant  on  all  flight  critical  sensors  and  control  effectors. 

In  order  to  support  design,  analysis,  and  simulation  of  reconfigurable  control  laws  for  the 
TAFA  aircraft,  a  math  model  was  created  using  the  structure  shown  in  Figure  A.2.  The  structure 
of  this  model  is  used  for  both  linear  analysis  and  nonlinear  simulation.  All  components  are  used 
for  linear  analysis.  Some  of  the  high  frequency  effects  are  excluded  from  the  nonlinear 
simulation  to  reduce  throughput  requirements. 
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Figure  A.2:  Block  Diagram  of  TAFA  Analysis  /  Simulation  Model 

The  components  of  Figure  A.2,  with  the  exception  of  the  flight  control  system,  describe  the 
TAFA  airframe  and  hardware  subsystems.  Descriptions  of  these  models  are  provided  in  the 
following  sections. 

A.lMass  Properties  Models 

The  TAFA  mass  properties  model  generates  the  vehicle  weight,  inertia  matrix,  center  of 
gravity,  and  inverse  inertia  nrorix  for  use  by  the  aerodynamics,  propulsion,  and  equations  of 
motion  modules.  A  static  model  (i.e.  mass  properties  do  not  change  with  fuel  bum)  is  used. 

The  model  supports  a  variety  of  vehicle  configurations  representing  different  stores  and  fuel 
loadings.  These  are.useful  in  evaluating  a  flight  control  design’s  sensitivity  to  knowledge  of  the 
vehicle’s  mass  properties.  The  mass  property  models  for  the  TAFA  were  taken  from  the  Boeing 
ASTOVL  program,  where  a  similar  configuration  was  explored.  The  TAFA  mass  property 
models  are  provided  in  Table  A.l. 
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Configuration 

Weight 

Ob) 

1\\ 

Ivy  I// 

(slu«-ft2) 

l\/ 

1  Sc  <4 

1)1. eg  Wl.cg 
(in) 

Clean,  2  AIM-9  &  2 

28064 

14578 

93604 

103855 

559 

373.5 

0 

100 

AIM- 120  internal  cany 

VL  with  6  AIM-120 

30494 

22329 

97090 

114520 

-234 

379 

0 

100 

Take  Off,  2  AIM-9, 

39684 

28154 

110616 

134103 

994 

386 

0 

100 

2  AIM- 120 

Take  Off,  2  GBU-24 

43714 

45624 

115963 

115963 

-1502 

393 

0 

100 

VL,  1  HARM 

28448 

19195 

95202 

109372 

90 

376 

-3.23 

100 

CTOL,  1  1200  lb  store 

28664 

19785 

95282 

110000 

36 

377 

-4.95 

100 

Combat  Wt,  2  AIM-9, 

34569 

19245 

99692 

114355 

682 

379 

0 

100 

2  AIM- 120,  60%  fuel 

Table  A.1:  TAFA  Mass  Properties  Models 


The  center  of  gravity  parameters  are  given  in  terms  of  a  fuselage  station  (FS),  butt  line 
station  (BL)  and  waterline  station  (WL).  These  quantities  are  measured  positive  aft  from  the 
nose,  out  the  right  wing,  and  up,  respectively. 


A.2  Aerodynamics  Models 

The  TAFA  aerodynamics  model  computes  the  body  axis  aerodynamic  forces  and  moments 
acting  on  the  vehicle.  Non-dimensional  aerodynamic  coefficients  are  computed  and 
dimensionalized  based  on  the  vehicle’s  flight  condition  and  geometry  to  create  the  forces  and 
moments.  The  non-dimensional  coefficients  are  a  function  of  the  vehicle’s  rigid  body  states, 
atmospheric  conditions,  and  control  surface  positions.  The  moment  coefficcucs  are  translated  to 
the  vehicle’s  center  of  gravity.  Reductions  in  control  surface  effectiveness  due  to  structural 
flexibility  are  also  modeled. 

The  TAFA  aerodynamics  model  is  based  on  Boeing  wind  tunnel  testing  of  a  tailless 
ASTOVL  configuration.  Control  increments  for  the  leading  edge  passive  porosity  and  forebody 
blowing  were  estimated  from  test  data  derived  on  other  Boeing  programs.  The  resulting 
database  provides  a  six  degree  of  freedom  aerodynamics  model  that  covers  the  complete  range  of 
flight  operations  from  Mach  0  to  2.5,  altitude  from  0  to  80,000  ft,  angles  of  attack  from  -2  to  48 
deg,  and  sideslip  angles  from  -20  to  20  deg.  The  nonlinear  database  includes  static,  dynamic 
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coupling  and  control  surface  effectiveness  terms.  The  aerodynamic  effects  of  external  stores  are 
not  modeled. 

The  baseline  aerodynamics  model  computes  the  non-dimensional  force  and  moment 
coefficients  at  the  reference  center  of  gravity  using  Equation  (A.1). 

C( )  =  C(\  (^>a»P)+X^(  )6i  r,  +X^(  P™  (A-l) 

i  *J 

The  form  of  this  equation  is  the  same  for  the  force  and  moment  coefficients,  and  thus  the  (  ) 
subscript  denotes  any  selected  force  or  moment  coefficient:  x,y,z,  for  the  axial,  lateral,  and 
vertical  forces,  ^,m,n  for  the  roll,  pitch,  and  yaw  moments.  The  terms  in  Eq.  (A.1)  are  looked  up 

from  tables  using  multi-dimensional  linear  interpolation. 

In  Equation  (A.1),  C(  ^  (M, a,  P  )  is  the  controls  neutral  term  which  represents  the  static 

forces  and  moments  acting  on  the  vehicle  when  all  of  the  control  effectors  are  at  their  zero 
positions.  This  term  is  a  function  of  the  vehicle’s  mach  M ,  angle  of  attack  a  and  sideslip  p . 

The  Xc(  )b  )kf/r,  term  represents  the  aerodynamic  forces  and  moments  due  to 

l 

deflections  of  the  individual  control  effectors.  The  control  increments  are  modeled  in  terms  of 
individual  surface  deflections  (e.g.  left  and  right)  rather  than  collective  and  differential 
deflections  to  support  subsequent  integration  of  faihre  and  damage  models.  The  control 
increments  are  scaled  by  a  factor  Kf/Ri  which  represents  a  loss  in  control  effectiveness  due  to 

structural  flexibility  of  the  surface  and  supporting  structure.  This  effect  will  be  discussed  in 
detail  in  the  following  section. 

An  additional  term,  ^C(  ^  (A/,a,6,,5;),  models  the  forces  and  moments  induced  by  the 

i.j 

interaction  between  control  effectors.  The  deflection  of  one  surface  may  alter  the  airflow  over 
another  surface,  thus  altering  its  effectiveness.  This  effect  is  modeled  by  the  controls  interaction 
term. 

The  final  term  in  the  aerodynamic  coefficient  build  up  is  the  dynamic  term,  which  consists  of 
the  aerodynamic  damping  derivatives  (e.g.  Cn  )  scaled  by  their  respective  rotational  rates  (e.g. 
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r ).  These  terms  are  scaled  by  »  where  /  is  the  reference  length  and  V  is  the  vehicle  s 

airspeed.  The  reference  length  is  the  mean  aerodynamic  chord  c  for  the  longitudinal  plane,  and 
the  wing  span  b  for  the  lateral  directional  plane. 

The  control  powers  of  the  various  effectors  are  compared  in  Figure  A.3  at  a  subsonic  flight 
condition  (Mach  0.6).  In  these  plots,  the  line  types  (solid,  dashed,  dotted)  on  the  plots 
correspond  to  the  order  that  the  effectors  are  listed  in  tbx  y-axis  label  of  the  plot. 

The  first  row  of  plots  compares  the  nose  up  and  no.se  down  pitch  control  power  of  the  trailing 
edge  flaps  and  canards  as  a  function  of  angle  of  attack.  The  pitching  moment  increment 
produced  by  each  control  effector  is  plotted  for  maximum  surface  deflections.  The  trailing  edge 
flaps  provide  more  nose  up  control  power  across  the  angle  of  attack  regime  due  to  the  small 
positive  deflection  limit  on  the  canard.  On  the  other  hand,  the  canard  provides  more  nose  down 
control  power  than  the  trailing  edge  flap  over  most  of  the  angle  of  attack  envelope. 

The  second  row  of  plots  compares  the  roll  control  power  of  the  ailerons,  trailing  edge  flaps, 
and  leading  edge  passive  porosity.  The  ailerons  and  trailing  edge  flaps  are  used  differentially  in 
this  case.  The  first  plot  compares  the  rolling  moment  coefficient  increment  generated  by 
maximum  deflection  of  each  effector  as  a  function  of  angle  of  attack.  The  trailing  edge  flaps 
have  the  most  control  power;  however,  this  analysis  assumes  that  all  of  the  deflection  capability 
is  available  for  roll.  In  reality,  the  differential  trailing  edge  flap  deflection  capability  may  be 
limited  based  on  the  collective  deflection  utilized  for  pitch  (i.e.  the  axis  prioritization  assigned  to 
the  control  effector).  The  aileron  provides  more  roll  control  power  than  the  leading  edge  passive 
porosity.  At  angles  of  attack  above  10  deg,  the  leading  edge  passive  porosity  provides  up  to  half 
of  the  aileron  control  power,  thus  making  it  a  viable  low  rate  roll  device.  The  second  plot 
compares  the  linearity  of  the  control  powers  as  a  function  of  deflection  (percentage  of  maximum) 
at  a  fixed  angle  of  attack.  The  aileron  and  trailing  edge  flap  control  powers  are  piecewise  linear 
with  a  reduced  control  derivative  at  higher  deflections.  The  leading  edge  passive  porosity 
control  power  is  linear  with  deflection. 

The  third  row  of  Figure  A.3  compares  the  yaw  control  power  of  the  aft  body  split  flaps, 
canards,  and  forebody  blowing.  The  canards  are  used  differentially  in  this  case. .  The  first  plot 
compares  the  yawing  moment  coefficient  increment  generated  by  maximum  deflection  of  each 
effector  as  a  function  of  angle  of  attack.  The  aft  body  split  flaps  are  most  effective  at  low  angles 
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of  attack,  while  the  forebody  blowing  is  most  effective  at  high  angles  of  attack.  The  second  plot 
compares  the  linearity  of  the  control  powers  as  a  function  of  deflection  (percentage  of  maximum) 
at  a  fixed  angle  of  attack.  The  aft  body  split  flap  and  forebody  blowing  control  powers  are  linear 
with  deflection,  while  the  differential  canard  control  derivative  varies  with  deflection.  The 
largest  control  derivative  occurs  in  the  low  to  moderate  deflection  region. 
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Nose  Up  Control,  Max  Deflections,  M=0.6 


Nose  Down  Control,  Max  Deflections,  M=0.6 


Roll  Control  for  Max  Deflections,  M=0.6 


Roll  Control  vs  Deflection,  M=0.6,  AoA=20  deg 


Yaw  Control  for  Max  Deflections,  M=0.6 


Yaw  Control  vs  Deflection,  M=0.6,  AoA=20  deg 


Figure  A 3:  TAFA  Aerodynamic  Control  Increments 
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Once  the  aerodynamic  force  and  moment  coefficients  are  computed  at  the  reference  center  of 
gravity,  the  aerodynamics  moments  are  translated  to  the  actual  center  of  gravity  of  the  vehicle. 
This  consists  of  calculating  the  incremental  moments  induced  by  the  aerodynamic  forces  acting 
through  a  point  displaced  from  the  center  of  gravity.  The  computation  is  shown  in  Equation 
(A.2). 
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(A.2) 


The  aerodynamic  force  and  moment  coefficients  are  then  dimensionalized  to  form  forces  and 
moments  using  the  vehicle  geometry  and  dynamic  pressure.  The  force  and  moment  equations 
are  summarized  in  Equation  (A.3). 

Fx=qSCx  L  =  qSbC, 

Fy  =  qSCy  M  =  qScCm  (A.3) 

Fz=qSCz  N  =  qSbCn 

The  forces  (Fx,  Fy,  Fz)  are  in  lbs,  and  the  moments  (L,  M,  N)  are  in  ft-lbs. 

The  airframe  geometry  constants  for  the  TAFA  aircraft  are  provided  in  Table  A.2. 


TAFA  Airframe  Geometry  Constants 

FSref  =  382  in 

S  =  590  ft2 

BLref=  0  in 

b  =  37.98  ft 

WLref=  100  in 

c  =  19.415  ft 

Table  A.2:  TAFA  Airframe  Geometry  Constants 


Nonlinear  Control 
Of  Fighter  Aircraft 

A.2.1  Aerodynamic  Flexibility  Effects 

The  aerodynamic  flexibility  effects  model  a  reduction  in  control  effectiveness  due  to 
structural  flexibility  of  the  control  surface  and  supporting  structure.  Derivation  of  these  effects 
requires  a  detailed  analysis  of  the  aerodynamics  loads  and  the  structural  characteristics  of  the 
vehicle.  This  data  was  not  available  for  the  TAFA  aircraft  and  thus  the  aerodynamic  flexibility 
effects  from  the  F- 15  ACTIVE  aircraft  were  used. 

The  aerodynamic  flexibility  model  reduces  the  control  effectiveness  as  a  result  of  wing  / 
surface  flexure  under  the  aerodynamic  loads.  The  actuator  takes  the  surface  to  its  commanded 
position;  however,  the  “effective”  position  (relative  to  the  air  stream)  is  less  due  to  flexing  of  the 
wing  and/or  surface.  This  effect  is  modeled  as  a  scale  factor  on  the  control  surface  aerodynamic 
increments. 

For  the  TAFA  planform,  the  aerodynamic  flexibility  effects  are  most  prominent  on  the 
trailing  edge  flap  and  aileron  effectors.  The  aft  body  split  flaps  consist  of  two  panels  (on  each 
side  of  the  aircraft),  one  of  which  opens  above  the  wing  and  the  other  below  the  wing.  The 
panels  tend  to  produce  torsional  moments  on  the  wing  in  opposite  directions.  These  torsional 
moments  cancel  each  other  to  nullify  wing  flexure  due  to  their  aerodynamic  loads  on  the  aft  body 
split  flap  panels.  Aerodynamic  flexibility  of  the  canards  is  not  expected  to  be  "significant  due  to 
their  low  aspect  ratio  and  moderate  leading  edge  sweep.  These  factors  tend  to  produce  a  small 
angular  deviation  of  the  elastic  axis  of  the  surface  from  the  vehicle  centerline,  thus  minimizing 
changes  in  the  local  angle  of  attack  of  the  surface  due  to  flexure.  Aerodynamic  flexibility  does 
not  affect  the  leading  edge  passive  porosity  or  forebody  blowing  since  these  effectors  since  they 
do  not  utilize  hinged  surfaces  operating  in  the  air  stream. 

As  a  result,  the  aerodynamic  flexibility  effects  are  applied  to  the  trailing  edge  flaps  and 
ailerons  only.  The  control  power  reduction  is  reflected  in  all  axes,  and  is  characterized  by  the 
flex-to-rigid  scale  factor  (KF/r)  shown  in  Figure  A.4.  The  scale  factor  is  commonly  a  function  of 
mach  and  altitude  to  reflect  movement  in  the  surface  center  of  pressure  with  mach,  as  well  as 
variations  in  loads  with  changes  in  dynamic  pressure. 
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Figure  A.4:  TAFA  Aerodynamic  Flexibility  Effect  Model 


A3  Propulsion  Model 

The  TAFA  propulsion  system  model  computes  the  body  axis  forces  and  moments  induced  by 
the  engine’s  gross  thrust,  the  ram  drag  of  the  engine  inlet,  and  the  nozzle  thrust  vectoring  angles. 
The  model  is  static  and  is  a  function  of  the  power  lever  angle  (PLA),  the  aircraft  rigid  body  states 
(rotational  rate,  angle  of  attack,  sideslip,  and  velocity),  and  the  aircraft  mass  properties  and 
geometry  data.  Nozzle  dynamics  are  modeled  in  the  actuation  system. 

The  propulsion  model  is  based  on  tabulated  steady-state  performance  data  of  a  moderate 
bypass  ratio  turbofan  engine.  The  gross  thrust  and  ram  drag  data  is  interpolated  from  these 
tables  based  on  the  aircraft’s  Mach  number,  altitude,  and  PLA.  The  propulsion  data  is  shown  in 
Figure  A.5  for  reference. 


Gross  Thrust  (lb)  Gross  Thrust  (lb)  Gross  Thrust  (lb) 
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Figure  A-5:  TAFA  Propulsion  System  Data 
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The  engine  gross  thrust  (FG)  data  is  resolved  to  the  engine  axes  by  transforming  it  through 
the  thrust  vectoring  angles.  The  axi-symmetric  nozzle  deflection  is  modeled  via  pitch  (5-p/p) 

and  yaw  (8xvy)  vectoring  angles.  Sequential  orthogonal  transformations  through  the  yaw  and 
pitch  vectoring  angles  resolve  forces  from  engine  to  nozzle  axes,  as  illustrated  in  Figure  A.6. 


Figure  A.6:  Relationship  Between  Engine  and  Nozzle  Axes 


The  propulsive  forces,  in  engine  axes,  are  defined  in  terms  of  the  engine  gross  thrust  and 
vectoring  angles  as 

FXe  =  Fg  TVp  )  TVy  ) 

FYe  =  Fg  cos^yvjjsin^xvy)  .  (A.4) 

FZe  =  fg  sin^yvp) 

The  body  axis  propulsive  forces  and  moments  are  formed  from  gross  thrust  and  inlet  (ram) 
drag  contributions.  The  gross  thrust  forces  are  transformed  from  engine  axes  to  body  axes  via 
rotations  through  the  engine  installation  angles.  These  angles  represent  the  orientation  of  the 
engine  with  respect  to  the  vehicle  body  axes.  The  transformation  from  body  axes  to  engine  axes 
is  accomplished  via  sequential  rotations  through  the  engine  yaw  ( \\r  E )  and  pitch  (  0  E ) 
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installation  angles.  The  roll  orientation  of  the  engine  is  assumed  to  be  coincident  with  the  body. 
The  body  axis  forces  due  to  gross  thrust  are  defined  as 

Fxg1  MQeWve)  -sin(¥E)  sin(0E)cos(\j/E)TFXE 

FYg  =  cos(0E)sin(\|/E)  cos(\yE)  sin(eE)sin(\|/E)  FYe  .  (A.5) 

fZg  [  -  sin(eE)  o  cos(eE)  J[fZe_ 

Moments  due  to  gross  thrust  are  computed  based  on  the  gross  thrust  forces  and  the  offset 
between  the  thrust  centerline  and  the  vehicle  center  of  gravity.  This  yields 

MLo  -(FSNoz-FSCg)/12  FXo 

MMfl  =  ?xFg  =  (BLN0Z-Blcg)/12  x  FY(,  .  (A.6) 

_MNg  J  L-(WLnoz  -■ WLcg)/12j  |_FZq  _ 

The  ram  drag  forces  and  moments  are  comprised  of  a  static  and  a  dynamic  contribution.  The 
static  forces  and  moments  are  dependent  on  the  vehicle’s  angle  of  attack  and  sideslip,  while  the 
dynamic  forces  and  moments  depend  on  the  vehicle’s  rotational  rates. 

The  ram  drag  static  forces  are  computed  from  the  inlet  ram  drag  (Dram)  and  the  vehicle’s 
angle  of  attack  and  sideslip  as 


"“Dram  *  cos(a)  *  cos(p)" 
-Dram*sin(|3) 
.“Dram  *  sin(cx)  *  c°s((3)_ 


(A.7) 


The  static  ram  drag  moments  are  based  on  the  ram  drag  forces  and  the  inlet  location  relative  to 
the  vehicle’s  center  of  gravity,  and  are  computed  as 
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The  dynamic  ram  drag  forces  are  based  on  the  vehicle’s  ram  drag,  rotational  rates  (p,  q,  r) 
and  airspeed  (V),  and  are  computed  as 
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(A.9) 


The  dynamic  ram  drag  moments  are  computed  from  the  dynamic  ram  drag  forces  based  on  the 
location  of  the  inlet  relative  to  the  vehicle’s  center  of  gravity.  The  computations  are  summarized 
below. 
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The  total  propulsive  forces  and  moments  are  formed  by  summing  the  gross  thrust,  static  ram 
drag,  and  dynamic  ram  drag  contributions.  The  TAFA  propulsion  mode  constants  are  provided 
in  Table  A.3. 


TAFA  Propulsion  Model  Constants 

FS„oZ  =  587.61  in 

FSram  =  208.87  in 

BLnoz  =  0  in 

BLram  =  0in 

WLnoz =  1 00  in 

WLram  =  83.23  in 

FSini  =  208.87  in 

6e  =  0  deg 
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BLini  =  0  in 

\|/e  =  0  deg 

WLjfli  =  83.23  in 

Table  A.3:  TAFA  Propulsion  Model  Constants 


A.4  Rigid  Body  Dynamics 

The  TAFA  Rigid  Body  Dynamics  model  implements  the  rigid  body  six  degree  of  freedom 
equations  of  motion.  The  translational  dynamics  are  propagated  in  an  inertial  reference  frame, 
while  the  rotational  dynamics  are  propagated  via  quaternions.  Additional  dynamic  parameters 
are  computed  and  output  from  the  model.  The  model  supports  flat  earth  equations  of  motion. 
The  inputs  to  the  model  are  the  total  body  axis  forces  and  moments  acting  on  the  vehicle,  as  well 
as  the  vehicle  mass  properties. 

The  TAFA  Rigid  Body  Dynamics  model  first  computes  the  translational  (linear)  and 
rotational  accelerations  acting  on  the  vehicle.  To  compute  the  linear  acceleration,  the 
acceleration  due  to  specific  forces  is  first  computed  as 

fB/m  (A.  11) 

where  m  is  the  vehicle  mass.  This  is  the  total  acceleration  due  to  external  forces,  excluding 
gravity.  Note  that  the  specific  force  vector  is  expressed  in  the  body  frame.  The  specific  force  is 
transformed  to  the  inertial  frame  and  added  to  gravity  to  form  the  total  linear  acceleration  of  the 
vehicle,  expressed  in  the  inertial  frame. 

a'  =  g'+C;(f5/m)  (A.  12) 


The  body-to-inertial  transformation  matrix  C'B  is  computed  from  the  body-to-inertial  attitude 


quaternion  q^  =  (q0, 


2  ,  2  2  2 
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(A.  13) 


Since  the  body-to-inertial  transformation  is  available  here,  the  gravity  vector  is  also  transformed 
to  the  body  frame  for  use  elsewhere. 
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g»=C?g'=(c;)V  (A>4) 

The  angular  acceleration  of  the  vehicle,  expressed  in  the  body  frame,  is  computed  as 


ci>*  = 


=  r1(m,-©ixl©i') 


(A.  15) 


where  I  is  the  inertia  matrix. 

The  vehicle  linear  and  angular  accelerations  are  then  integrated  to  produce  velocity,  position, 
angular  rate,  and  attitude.  The  linear  inertial  acceleration  vector  is  integrated  to  produce  a  linear 
inertial  velocity  vector. 

v'  =  j  a'dt  (A.16) 

The  linear  inertial  velocity  vector  is  integrated  to  produce  an  inertial  position  vector. 


r/  =  J \'dt  (A.  17) 

The  angular  acceleration  vector  is  integrated  to  produce  an  angular  velocity  vector. 

>" 

\  =  \&Bdt  (A.  18) 


t o*  = 


Using  the  body  attitude  quaternion,  the  angular  velocity  vector  is  transformed  to  form  an  attitude 
quaternion  rate 
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(A.  19) 


which  is  integrated  to  produce  a  body  attitude  quaternion 


q* 


(A.20) 


In  all  cases,  the  preceding  vector  integrals  are  reducible  to  independent  scalar  integrals  on  die 
components  of  the  vectors.  The  scalar  integrals  are  evaluated  numerically  using  the  integration 
method  and  time  step  selected  by  the  user. 

The  TAFA  Rigid  Body  Dynamics  model  then  converts  the  vehicle  state  vector  elements  from 
the  inertial  frame  to  a  geographic  (local  level)  frame  according  to  the  Earth  model.  The  Flat 
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Earth  model  assumes  that  the  inertial  coordinate  frame  used  in  propagating  the  state  vector  is  the 
geographic,  North-East-Down  coordinate  frame.  It  simply  passes  the  position  and  velocity 
vectors  through,  renaming  their  elements  to  North,  East,  Down,  and  Vn,  Ve,  Vd  respectively. 
The  body-to-inertial  attitude  quaternion  is  also  passed  through  and  renamed  to  the  body-to- 
geographic  attitude  quaternion.  The  position  vector  is  transformed  to  Latitude,  Longitude, 
Altitude  (Latitude  and  Longitude  are  computed  from  the  North  and  East  positions  divided  by  the 
earth  radius,  while  Altitude  is  computed  as  -Down).  Finally,  the  gravity  vector  for  the  Flat  Earth 
model  is  computed,  which  is  just  a  constant  [0, 0,  32.174]  ft/s/s. 

Once  the  dynamics  propagation  is  completed,  additional  useful  quantities  are  computed 
which  may  be  required  as  flight  control  feedback  or  inputs  to  sensor  models.  The  North-East- 
Down  components  of  the  vehicle  velocity  are  used  to  compute  the  vehicle's  speed  (velocity 
magnitude)  and  the  flight  path  angles.  The  speed  is  computed  as 

Velocity _frs  =  \  |vf°|l  =  +Vl  +^j)  (A.21) 


The  vertical  flight  path  angle  is  computed  as 

(  -  \ 


y  =  sin 


-i 


~Fn 


(  Velocity _fps 


(A.22) 


while  the  heading  (horizontal  flight  path)  angle  is  computed  as 


X  =  tan 


-i 


\V»/ 


(A.23) 


The  wind  velocity  v^£°  at  the  vehicle's  location,  which  is  the  velocity  of  the  local  air  mass, 
is  used  to  compute  the  vehicle  velocity  with  respect  to  the  air. 

„NED  _  -.NED  -.NED 

VSWr  =  Vfi  ~vw  .  (A.24) 

This  air-relative  velocity  is  transformed  to  the  body  frame  to  form  the  body-frame  air-relative 
velocity, 

r* 

U_wrt_air 

(A.25) 


^  B.air 


V_wrt_air 
W  wrt  air 


which  is  subsequently  used  to  compute  aerodynamic  angles  and  rates. 
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The  geographic-to-body  transformation  matrix  C*G  is  computed  from  the  body-to-geographic 
transformation  quaternion  q°B  =(q0,qltq2,q3)-  A  matrix  transpose  is  included  in  the  calculation 
to  produce  a  matrix  which  will  transform  in  the  opposite  sense  as  the  quaternion. 

~q\  +  q]  -  q\  -  q]  2 {q#2  +  q0q2 )  2(q]q3  -  q0q2 ) 


CB  = 

'-'G 


(A.26) 


2  (qxq2 -q,qi)  q\~q\  +  *2  “ft*  +«ofi) 

2 (<7^3  +  )  2(9^3  -  )  q\  -  q]  ~  q\  +  q\ 

The  geographic  frame  is  the  same  as  the  North-East-Down  frame,  so  the  geographic-to-body 

transformation  matrix  CBG  can  be  used  as  a  North-East-Down-to-body  transformation  matrix 

Cned  •  Tbe  matrix  is  used  to  transform  the  body  velocity  and  the  body  air-relative  velocity  to  the 


body  frame. 
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The  rate  of  change  of  the  body  components  of  velocity  are  computed  using  the  Coriolis  equation, 
v*  =  (dv /  dt))B  -ffl'xvj 


U 

V 

W 


=  aB  ~(bB  x  \bb 


(A.28) 


where 

aB  =  fB/m+gB.  (A-29) 

The  Euler  angles  describing  the  orientation  of  the  body  relative  to  the  local  level,  North-East- 
Down  frame  are  computed  from  elements  of  the  geographic-to-body  transformation  matrix. 
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Roll  =  o  =  ATANl{CBG  (2,3), C*  (3,3)) 

Pitch  =  0  =  sin'1  (-  CBc  (13))  (A.30) 

Yaw='¥  =  ATAN2(CBG(1,2),CBG(.U)) 

The  ATAN2  function  is  the  four-quadrant  arctangent  function.  To  prevent  an  argument  error 
due  to  numerical  round  off  errors,  the  argument  of  the  arcsine  function  is  limited  to  the  range  ±1 
inclusive 

The  vehicle's  air-relative  velocity  and  its  acceleration  and  angular  rate  are  used  to  compute 
the  true  airspeed,  aerodynamic  angles,  and  aerodynamic  angle  rates.  Two  of  the  aerodynamic 
angles  (angle  of  attack  and  sideslip)  are  then  used  to  transform  the  body  angular  velocity  vector 
to  the  wind  axes  coordinate  frame.  Because  the  wind  frame  is  a  rotating  frame,  this  is  not  the 
body  angular  velocity  with  respect  to  the  wind  frame,  but  merely  the  projection  of  the  body 
angular  velocity  into  the  wind  frame. 


V 

cos  p  cosa 

sinP 

cos  p  sina 

>" 

a 

-B 

=  cB©  = 

-  sin  P  cosa 

cosP 

-  sin  P  sina 

Q 

_RW_ 

-  sin  a 

0 

cosa 

R 

The  true  airspeed  is  computed  as 

^tas  =  Airspeed — fps  —  ||vjiaj;.|  .  (A. 32) 

The  aerodynamic  angles  are  computed  as  follows. 


jt  =  ATAN2(W,U) 

P  =  sin'1  (v  /  VTAS  ) 

Pfl  =  ATAN2(V,U) 
ar  =cos  ~'(U/VTAS) 
aT  =  cos~'({/  /  VTAS  ) 


Angle  of  attack 
Sideslip  angle 
"Body"  Sideslip  angle 
Total  angle  of  attack 
Aerodynamic  roll  angle 


To  prevent  an  argument  error  due  to  numerical  round  off  errors,  the  argument  of  the  inverse  sine 
and  inverse  cosine  functions  are  limited  to  the  range  ±1  inclusive.  If  Airspeed  is  zero,  then 
sideslip  and  total  angle  of  attack  are  set  to  zero. 
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The  aerodynamic  angle  rates  are  undefined  if  the  airspeed  is  zero.  Furthermore,  the  angle  of 
attack  rate  and  sideslip  rate  are  undefined  if  the  sideslip  angle  is  ninety  degrees  (causing 
U=W=0),  and  the  total  angle  of  attack  rate  is  undefined  if  the  total  angle  of  attack  is  zero  (so 
V  =  W  =  0 ).  The  rates  are  set  to  zero  in  these  circumstances.  However,  rather  than  performing 
an  explicit  check  that  values  are  equal  to  zero  and  setting  the  results  to  zero,  a  small  value  is 
added  to  the  airspeed  and  to  W2 .  Both  of  these  quantities  are  greater  than  or  equal  to  zero; 
adding  a  small  value  to  each  makes  them  strictly  greater  than  zero. 

vTASl  =  vTAS  +  ht6 
w2  =  w2+  nr6 

Using  these  perturbed  quantities  in  the  denominators  of  the  equations  for  aerodynamic  angle 
rates  prevents  division  by  zero  errors  and  gives  results  of  zero  in  the  cases  described  above.  The 
error  introduced  is  negligible. 

The  aerodynamic  angle  rates  are  computed  as 

.  UW-WU 


a  = 


U2+W2 


p  = 


ar  = 


VusV-IZm 


^TASl 


(A.33) 

(A.34) 

(A.35) 


when 


yTAs  = 


o v,v,w)*{u,v,w ) 


TAS 1 


and 


(0,  V,  W)  =  (f  /  m  +  g)*  -  (5  *  x  vL,  (assumes  wind  is  not  time-varying) 


(A.36) 


(A.37) 


Finally,  load  factors  are  computed  from  the  specific  forces  acting  on  the  vehicle.  The  body- 
frame  Normal  and  Side  load  factors  are  computed  from  the  body-frame  Z  and  Y  components  of 
specific  force  respectively. 
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NZCG_g_  plus_  up  =  - 


f B(z)/ m 
32J74 


(A.38) 


NYCG_g  = 


32.174 


(A.39) 


The  Normal  Load  Factor  is  computed  from  the  stability  frame  Z  component  of  specific  force, 
i.e.,  the  acceleration  due  to  Lift  force.  The  Lift  force  is  found  by  transforming  body-frame  X  and 
Z  components  of  specific  force  by  angle  of  attack. 


a  s{z)  =  L/m  =  /  m)sina  +  (f  B  (z)  /  ro)cosa 

as(z) 

“J-Tit  . 


(A.40) 

(A.41) 


A.5  Flexible  Vehicle  (Structural)  Dynamics 

The  flexible  vehicle  dynamics  model  captures  the  effects  of  the  vehicle’s  structural  flexibility 
on  the  sensor  measurements.  These  effects  are  dependent  on  the  location  of  the  sensors  and  the 
shapes  of  the  various  structural  modes.  The  influence  of  the  structural  modes  on  the  sensor 
measurements  must  be  attenuated  through  the  use  of  aeroservoelastic  (ASE)  compensation  filters 
prior  to  use  in  the  control  law,  since  insufficient  attenuation  of  these  effects  can  lead  to 
instability. 

Linear  analysis  is  commonly  used  to  evaluate  the  impact  of  the  structural  dynamics  on  the 
control  system  performance  and  to  determine  the  ASE  filtering  requirements.  Linear  models  of 
the  ASE  dynamics  are  generated  at  various  flight  conditions  using  an  integrated  high  fidelity 
model  of  the  vehicle’s  aerodynamic  and  structural  characteristics.  The  model  determines  the 
excitation  of  the  structural  dynamics  caused  by  the  aerodynamic  and  inertial  loads  acting  on  the 
vehicle,  and  the  subsequent  impact  on  sensor  measurements.  The  model  also  captures  the  effect 
of  the  vehicle’s  flexure  on  the  rigid  body  dynamics.  The  resulting  linear  ASE  dynamics  model  is 
then  combined  with  linear  models  of  the  rigid  body  dynamics,  control  system,  and  avionics 
subsystems  to  design  ASE  compensation  filters  and  evaluate  system  performance. 

The  6DOF  model  captures  the  effects  of  the  structural  flexibility  on  the  sensor 
measurements,  but  neglects  the  influence  of  the  structural  flexibility  on  the  rigid  body  dynamics 
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since  this  is  generally  a  secondary  effect  The  equations  used  to  model  the  structural  dynamics 
are  based  on  those  described  in  References  2  and  3. 

The  equations  governing  the  structural  mode  dynamics  have  the  form 

[s2  +2^,(S)ls+(iif%  =  —  X M5 j M5; )+ fa' (5; )m*j i R>  +a<(5yKy]s25^^]  (A-42) 

mi  j  L 


where 


£  -  Generalized  coordinate  of  the  i*  structural  mode 
£  (  -  Damping  of  the  i*  structural  mode 
(D(.  -  Natural  frequency  of  the  i*  structural  mode  (rad/sec) 
m,  -  Generalized  mass  of  the  i*  structural  mode  (slug) 

<f)/  (Sj )  -Deflection  of  the  i*  structural  mode  at  the  location  of  control  effector  8 }  (ft/ft) 

j  -Force  due  to  control  effector  5 y  which  acts  to  excite  the  structural  mode  of  interest. 

This  force  lies  in  the  plane  of  interest  for  bending  modes.  For  torsional 
modes,  the  direction  of  the  moment  producing  force  is  used. 

a.  ^  j  .  Slope  of  the  i*  structural  mode  at  the  location  of  control  effector  8  y  (rad/ft) 


mD 


-Mass  of  control  effector  8 ,  (slug) 


£  n  -Offset  between  the  center  of  gravity  and  the  hinge  line  of  control  effector  8 y  (ft) 

RJ 

I  -  Moment  of  ine~ua  of  control  effector  8y  about  its  hinge  line  (slug-ftJ) 


The  <j>,  (bj  term  on  the  right  hand  side  represents  the  excitation  of  the  i*  structural  mode 

by  the  direct  forces  applied  by  the  control  effectors.  The  remaining  terms  on  the  right  hand  side 
represent  the  effect  of  the  inertial  reaction  torque  induced  by  the  angular  acceleration  of  the 
control  effectors  on  the  excitation  of  the  structural  mode.  The  excitation  forces  are  summed  over 
the  control  effector  suite. 

The  effect  of  the  structural  modes  on  the  sensor  measurements  depends  on  the  location  of  the 
sensors  along  the  mode  shape.  The  effect  of  a  given  mode  is  proportional  to  the  slope  or 
deflection  of  the  mode  at  the  sensor  location.  The  overall  effect  of  the  vehicle’s  flexibility  on  the 
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sensor  measurements  is  the  sum  of  the  contributions  from  the  various  modes.  Sensed  quantities 
for  a  given  mode  are  related  to  the  mode  shape  as  follows 


a,  =a,(a£, 

Pz  =°/(P& 

Pn<lnri  =^i(Gyro%i 
hi  =  <j>, 

Ayi,Ah  =<t )t{Accel% 
4 =  G,{Gyro% 


(Angle  of  Attack) 
(Sideslip) 

(Rotational  Rates) 
(Altitude) 

(Linear  Acceleration) 
(Body  Attitudes) 


For  these  equations,  the  quantity  in  parenthesis  indicates  the  point  at  which  the  noted  mode  slope 
or  deflection  is  defined.  The  mode  slopes  and  deflections  are  defined  with  respect  to  the  axis  of 
interest. 

Since  detailed  structural  mode  characteristics  of  the  TAFA  aircraft  were  not  available,  a 
structural  model  was  created  by  combining  the  structural  mode  data  of  the  F- 15  ACTIVE  with 
the  TAFA  aerodynamics.  This  model  is  approximate  since  it  assumes  the  structural  mode 
frequencies  and  mode  shapes  for  the  TAFA  and  F- 15  ACTIVE  are  the  same,  but  retains  realism 
by  capturing  the  structural  mode  excitation  supplied  by  the  relative  forces  and  moments  exerted 
by  each  control  effector.  The  contribution  of  the  ASE  dynamics  to  the  sensor  measurements  was 
verified  to  be  comparable  to  that  of  the  F-15  ACTIVE. 

The  first  step  in  the  development  of  the  ASE  dynamics  model  was  to  identify  grid  points  on 
the  F-15  ACTIVE  structural  model  which  correspond  to  the  approximate  locadon  of  the  control 
effectors  and  sensors  on  the  TAFA  aircraft.  These  grid  points  are  required  to  define  the  mode 
slopes  and  deflections  used  in  the  ASE  dynamics  differential  equations.  Since  the  aircraft  are 
similar  in  size,  the  TAFA  control  effector  and  sensor  locations  were  assumed  to  correspond  to 
similarly  located  devices  on  the  F-15  ACTIVE  planform.  The  grid  points  used  are  provided  in 
Table  A.4. 

The  next  step  was  to  model  the  forces  and  inertial  reaction  torques  which  form  the  forcing 
terms  that  excite  the  structural  modes.  The  direct  force  terms  are  computed  directly  from  the 
TAFA  aerodynamics  and  propulsion  models,  and  are  simply  the  force  increments  due  to  the 
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various  control  effectors.  These  terms  come  directly  from  the  nonlinear  aerodynamic  and 
propulsion  models  and  thus  vary  with  flight  condition  and  maneuver  level. 


TAFA  Sensors  /  Control  Effectors 

F-15  ACTIVE  Sensors  /  Control  Effectors 

Trailing  Edge  Flap,  Aft  Body  Flap,  Aileron 

Stabiiaior 

Thrust  Vectoring  Nozzle 

Thrust  Vectoring  Nozzle 

Canard 

Canard 

Leading  Edge  Passive  Porosity 

Aileron 

Forebody  Blowing 

(Location  on  Forward  Forebody) 

Accelerometers 

Accelerometers 

Gyros 

Gyros 

AoA  /  Sideslip 

AoA  /  Sideslip 

Altimeter 

Accelerometer  Location 

Table  A.4:  Grid  Points  Used  in  Development  of  TAFA  Structural  Dynamics  Model 


The  inertial  reaction  torque  terms  are  a  function  of  the  geometry  of  the  control  effectors. 
Since  a  detailed  structural  model  of  the  TAFA  was  not  available,  the  control  effector  geometry 
data  was  estimated  as  described  below.  The  inertial  reaction  torque  only  exists  for  control 
effectors  which  rotate  a  significant  mass  about  a  pivot  point.  As  a  result,  the  leading  edge 
passive  porosity  and  forebody  blowing  do  not  produce  a  significant  inertial  reaction  torque.  The 
mass,  offset  between  the  effector  center  of  gravity  ard  pivot  axis,  and  the  moment  of  inertia  of 
the  effector  about  its  pivot  point  must  be  estimated. 

The  mass  of  the  aerodynamic  control  effectors  was  estimated  by  multiplying  the  surface  area 
of  the  control  effectors  by  a  density  (mass/surface  area).  The  density  was  computed  using  mass 
property  data  for  the  F/A-l  8C/D  control  surfaces,  and  assumed  to  be  the  same  for  the  TAFA 
aerodynamic  control  effectors.  The  density  of  the  F/A-l  8C/D  trailing  edge  flap  (TEF)  was  used 
for  the  subsequent  computations,  and  was  found  to  be: 


F/A- 1 8C/D  TEF  Weight  1 1 7.39  lb 

F/A- 1 8C/D  TEF  Surface  Area  47.45  ft2 

density  2.47  lb/ft2  =  0.07689  slug/ft2 
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This  led  to  the  mass  values  for  the  TAFA  aerodynamic  control  effectors  provided  in  Table  A.5. 
The  mass  of  the  thrust  vectoring  nozzle  was  based  on  mass  property  data  for  the  F- 15  ACTIVE, 
and  is  also  provided  in  the  table. 


Control  Effector 

Surface  Area  (ft1) 

Mass  (slug) 

Trailing  Edge  Flap 

25.2 

1.94 

Aft  Body  Split  Flap 

10.9 

0.84  | 

Aileron 

7.7 

0.59 

Canard 

20.1 

1.54 

Thrust  Vectoring  Nozzle 

- 

9.1 

Table  A.5:  Mass  Properties  of  TAFA  Control  Effectors 


The  offset  between  the  center  of  gravity  and  pivot  point  of  each  control  effector,  £ r.  ,  was 

estimated  from  the  geometry  of  the  control  effectors.  The  effectors  were  assumed  to  have 
uniform  density  and  thus  the  center  of  gravity  is  at  the  geometric  center.  The  distance  from  the 
geometric  center  to  the  pivot  point  was  thus  taken  as  the  required  offset.  The  values  are  provided 
in  Table  A.6. 


Control  Effector 

Offset  Between  eg 

and  Pivot  Point  (ft) 

Trailing  Edge  Flap 

1.68 

Aft  Body  Split  Flap 

1.28 

Aileron 

0.7 

Canard 

1.6 

Thrust  Vectoring  Nozzle 

1.0 

Table  A.6:  Offset  Between  eg  and  Pivot  Point  of  TAFA  Control  Effectors 
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The  final  quantity  required  for  the  inertial  reaction  torque  model  is  the  moment  of  inertia  of 
each  control  effector  about  its  pivot  axis.  This  data  was  estimated  from  the  mass  moments  of 
inertia  for  simple  geometrical  shapes  (e.g.  a  thin  rectangular  plate).  The  resulting  inertia  data  is 
provided  in  Table  A.7. 


Moment  of  Inertia 

Control  Effector 

(slug-ft2) 

Trailing  Edge  Flap 

36.38 

Aft  Body  Split  Flap 

5.2 

Aileron 

5.95 

Canard 

3.88 

Thrust  Vectoring  Nozzle 

75 

Table  A.7:  Moment  of  Inertia  of  TAFA  Control  Effectors 

Using  this  data,  the  structural  dynamics  were  integrated  into  the  TAFA  6DOF.  Longitudinal 
(symmetric)  and  lateral  directional  (anti-symmetric)  modes  were  implemented  separately.  The 
lateral  directional  model  includes  lateral  axis  bending  as  well  as  torsional  mode  effects.  These 
effects  combine  to  form  the  forcing  function  for  the  lateral  directional  modal  dynamics.  The 
ASE  model  in  the  TAFA  6DOF  is  configured  so  that  linear  models  can  be  extracted  at  selected 
operating  points. 

An  example  of  the  linear  ASE  dynamics  model  extracted  from  the  TAFA  6DOF  is  provided 
in  Figure  A.7.  This  figure  shows  the  magnitude  of  the  frequency  response  relating  the  ASE 
contribution  to  sensed  pitch  rate  and  normal  acceleration  to  the  pitch  axis  control  effector  inputs. 
This  data  is  for  the  low  altitude  ingress  /  egress  flight  condition  described  in  Reference  1.  This 
data  shows  high  frequency  peaks  in  the  frequency  response  curves  which  exceed  20  dB, 
indicating  the  need  for  ASE  filtering  in  the  control  law  design. 
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Figure  A.7:  Frequency  Respond  of  Longitudinal  Axis  Flexible  Vehicle  Dynamics  Model  for 

Low  Altitude  Ingress  /  Egress  Flight  Condition 


A.6  Actuator  Models 

The  actuator  models  describe  the  dynamic  response  characteristics,  position  limits,  and  rate 
limits  of  the  TAFA  control  effectors.  Both  high  and  low  fidelity  actuator  models  are  available. 
The  low  fidelity  model  assumes  position  and  rate  limit  capabilities  based  on  nominal  operating 
conditions.  The  high  fidelity  model  computes  the  position  and  rate  limit  capabilities  based  on 
the  surface  loads  (i.e.  hinge  moments).  The  hinge  moment  effects  reduce  the  position  and  rate 
capabilities  of  the  actuators  at  high  dynamic  pressure  conditions. 

The  sign  conventions  of  the  TAFA  control  effectors  are  described  in  Figure  A.8. 
Aerodynamic  effectors  and  thrust  vectoring  sign  conventions  are  described  in  terms  of  a  positive 
deflection  of  the  trailing  edge  (TF).  A  positive  command  to  the  forebody  blowing  effector  vents 
thrust  out  the  left  side  of  the  ah  craft  to  yield  a  positive  yawing  moment  Positive  commands  to 
the  leading  edge  passive  porosity  and  aft  body  split  flap  (left  or  right)  effectors  cause  them  to 
open  relative  to  their  neutral  position. 


Leading  Edge  Passive  Porosity 
(+  open) 

Canard 
(+  TE  Down) 


Axi-Symmetric  Pitch  /  Yaw  Thrust  Vectoring 

Aft  Body  Split  Flap  (+ open)  (+  TE  Down/+  TE  Left) 

Trailing  Edge  Flap  (+  TE  Down) 

Aileron  (+  TE  Down) 


"Forebody  Blowing 
(+  -*  +  Yaw) 

Figure  A.8:  TAFA  Control  Effector  Sign  Conventions 
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The  actuator  dynamics  of  the  various  control  effectors  are  modeled  using  first  or  second 
order  transfer  functions  with  position  and  rate  limits.  The  position  and  rate  limits  are  applied  to 
the  integrator  states  rather  than  the  integrator  outputs  to  prevent  windup  during  periods  of 
saturation.  The  first  order  actuator  dynamics  are  characterized  by  a  bandwidth  (toB),  while  die 

second  order  actuator  dynamics  are  characterized  by  a  damping  ratio  (£) and  natural  frequency 
( con).  The  linear  transfer  functions  of  these  models  are  provided  below. 

5  _  °>B  8  _  mn  (A.43) 

Sc  S  +  COb  5c  s2 +2£(0nS  +  C0n 

The  dynamic  parameters  used  to  characterize  the  TAFA  actuators  are  provided  in  Table  A.8. 


Surface 

Essmmm 

Canard 

LE  Passive  Porosity 
Trailing  Edge  Flap 
Aileron 

Aft  Body  Flap 
Vectoring  Nozzle 
Forebody  Blowing 

C0b  =  40  rad/s 
tOg  =  30  rad/s 

C0g  =  40  rad/s 
dig  =  40  rad/s 
€0g  =  40  rad/s 
£=0.8,  (D,  =  80  rad/s 

C0g  =  20  rad/s 

Table  A.8:  TAFA  Actuator  Dynamics 


The  low  fidelity  actuator  models  use  position  and  rate  limits  which  characterize  the 
capability  under  moderate  loads.  These  limits  are  provided  in  Table  A.9. 


Surface 

Position 

Limits 

Rate 

Limits 

Canard 

-80°/10° 

70°/s 

LE  Passive  Porosity 

0/1 

10  s-' 

Trailing  Edge  Flap 

-30°/45° 

90°/s 

Aileron 

-30730° 

1207s 

Aft  Body  Flap 

0790° 

1207s 

Vectoring  Nozzle 

-30730° 

60°/s 

Forebody  Blowing 

±0.006 

0.06  s’1 

Table  A.9:  TAFA  Low  Fidelity  Actuator  Model  Position  and  Rate  Limits 
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The  high  fidelity  actuator  models  reduce  the  position  and  rate  limit  capabilities  of  the 
actuator  due  to  the  aerodynamic  loads  acting  on  the  effector  at  the  current  operating  condition. 
Worst  case  models  are  formed  by  assuming  the  load  always  oppose  motion  of  the  surface.  The 
aerodynamic  loads  do  not  affect  the  leading  edge  passive  porosity,  forebody  blowing  and  thrust 
vectoring  effectors,  and  thus  their  position  and  rate  limits  are  unchanged  from  the  low  order 
model.  Position  and  rate  limit  capabilities  of  the  aerodynamic  control  effectors  are  characterized 
by  the  torque-speed  curve  shown  in  Figure  A.9. 


Rate  Limit  (8) 


Figure  A.9:  High  Fidelity  Actuator  Model  Torque-Speed  Curve 


The  torque-speed  curve  is  characterized  by  a  no  load  rate  limit  ( 5  n0  load )  ®  stall  torque 

(HMmax)-  The  stall  torque  is  the  maximum  load  which  the  actuator  can  oppose  to  move  the 
surface.  The  position  and  rate  limits  are  computed  from  the  torque-speed  curve  based  on  the 
current  hinge  moment  load  acting  on  the  surface.  The  torque-speed  is  characterized  by  the 
equation 

,2 


8  ] 

2+(  hm  ] 

no  load  > 

Ihm^J 

=  1  . 


(A.44) 


The  hinge  moment  of  each  surface  is  computed  based  on  a  dimensionalized  hinge  moment 
derivative  and  the  surface  deflection.  The  equation  is 

HM=q*Ssurf*f*CHM*|8|  (A.45) 


where 
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q  Dynamic  pressure 

5  surf  Surface  area  of  the  effector 

£  Offset  between  surface  center  of  pressure  and  hinge  line 
Chm  Hinge  moment  derivative  for  given  effector 

6  Deflection  of  given  effector. 


1  he  position  limits  for  a  given  effector  are  computed  as  the  minimum  of  the  deflections 
which  yield  the  associated  stall  torque  and  the  physical  deflection  limits  of  the  effector.  The 
equations  are  provided  below. 


HM 


max 


'““HM  q  *Ssurf  *Chm 


^  upper  limit  *n^^maxj^i  »  ^  max  physical  limit  ) 


'maxHM  ^maxptysicai 

^lower  limit  —  max(— 1 ^maxjjM  ’^““physical  limit ) 

The  corresponding  rate  limit  is  computed  from  the  equation  for  the  torque-speed  curve  as 


(A.46) 


(A.47) 


The  rate  limit  capability  is  limited  by  a  maximum  achievable  rate  (5^^)  which  characterizes 
the  actuator’s  capability  under  the  mechanical  load  imposed  by  the  fin  inertia,  gearing,  linkages, 
etc.  The  rate  limit  (5^^)  is  assumed  to  equal  the  rate  limit  capability  of  the  low  fidelity 
actuator  model  so  that  the  high  fidelity  model’s  rate  capability  will  not  exceed  that  of  its  low 
fidelity  counterpart. 

Data  was  obtained  from  Boeing’s  ASTOVL  program  indicating  stall  torque  and  no  load  rate 
limit  requirements  for  the  actuators.  These  values  are  assumed  for  the  TAFA  aircraft  and  are 
provided  in  Table  A.  10. 


No  Load  Rate 

Maximum  Rate 

Limit 

Capability 

Stall  Torque 

(^noloadl 

(Slim  it  1 

Surface 

(in-lb) 

(deg/s) 

(deg/s) 
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Trailing  Edge  Flap 

320000 

145 

90 

Aileron 

57000 

170 

120  1 

Canard 

290000 

105 

70 

Aft  Body  Flap 

35000 

180 

120 

Table  A.10:  TAFA  High  Fidelity  Actuator  Model  Parameters 

The  remaining  parameters  required  for  the  TAFA  high  fidelity  actuator  model  are  those  used 
to  compute  the  hinge  moments  of  the  various  control  effectors  as  a  function  of  dynamic  pressure 
and  surface  deflection.  A  detailed  hinge  moment  database  would  be  collected  by  instrumenting 
the  control  effectors  during  wind  tunnel  testing  and  measuring  the  shear  force,  root  bending 
moment  and  hinge  moment  as  a  function  of  flight  condition  and  surface  deflection.  This 
database  would  include  variations  with  mach  and  angle  of  attack  to  characterize  changes  in  the 
air  load  acting  on  a  given  surface  and  movement  of  the  center  of  pressure  of  the  surfaces.  A  high 
fidelity  hinge  moment  database  is  not  available  for  the  TAFA  aircraft  and  thus  a  simple  model 
will  be  used  to  capture  the  effects  of  hinge  moments  on  the  actuator  position  and  rate 
capabilities. 

Preliminary  data  obtained  from  the  ASTOVL  program  showed  that  the  trailing  edge  flaps 
and  ailerons  could  reach  their  full  deflection  at  dynamic  pressures  up  to  560  psf.  Comparable 
capability  was  assumed  for  the  canard  and  aft  body  flap,  and  used  to  compute  the  hinge  moment 
model  parameters.  The  surface  area  of  the  effectors  and  the  offset  between  the  center  of  pressure 
and  surface  hinge  line  were  estimated  from  the  TAFA  planform.  The  hinge  moment  derivative 
was  then  computed  assuming  that  the  maximum  surface  deflection  yielded  a  hinge  moment  equal 
to  the  stall  torque  when  the  dynamic  pressure  was  560  psf.  The  resulting  hinge  moment  model 
parameters  a^e  provided  in  Table  A.l  1 . 


Surface 

Surface  Area 

s  (ft2) 

Center  of 

Pressure  to 

Hinge  Line 

Offset 

/(in) 

Hinge  Moment 

Derivative 

Chm  (1/deg) 

Trailing  Edge  Flap 

25.2 

20.2 

0.0246 

Aileron 

7.7 

8.4 

0.0552 
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Canard 

202 

19.2 

0.0167 

Aft  Body  Flap 

10.9 

15.4 

0.0043 

Table  A.11:  TAFA  Hinge  Moment  Model  Parameters 


A.7  Sensor  Models 

The  TAFA  sensor  models  simulate  the  expected  signal  characteristics  yielded  by  the  sensors 
on  board  the  aircraft.  These  sensors  include  accelerometers,  gyros,  angle  of  attack  and  sideslip 
probes,  the  air  data  system,  and  the  navigation  /  attitude  reference  system.  All  sensed  signals  are 
computed  as  the  true  value  altered  by  a  time  delay,  measurement  bias,  and  zero  mean  white 
noise.  The  sensor  model  for  a  given  feedback  channel  is  shown  in  Figure  A.  10. 


True 

Signal 


Equations 
of  Motion 


Time 

Delay 


Bias  Gaussian 
White  Noise 


Sensed  Signal 


Figure  A.10:  TAFA  Sensor  Model 


Specifications  for  the  maximum  time  delay,  bias  level  and  noise  standard  deviation  for  each  type 
of  sensor  are  provided  in  Table  A.  12.  The  time  delay  parameter  includes  generation  of  the  signal 
as  well  as  the  time  between  signal  generation  and  utilization  by  the  flight  control  algorithms. 


Sensor  Type 

Time  Delay 

(msec) 

Bias 

Noise  Standard 

Deviation 

Accelerometers 

AX,  AY,  AZ  (ft/s") 

0 

0.1 

0.5 

Gyros 

p  (deg/s) 

0 

0.0 

0.32 

q,r  (deg/s) 

0 

0.0 

0.08 

AoA  /  Sideslip  Probes 
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Angle  of  Attack  (deg) 

0 

0.2 

0.2 

Sideslip  (deg) 

0 

0.1 

0.2 

Air  Data  System 

f  i 

Mach 

10 

0.01 

0.01 

Airspeed  (ft/s) 

10 

1.0 

10.0 

Dynamic  Pressure  (psf) 

10 

0.0 

10.0 

Navigation  /  Attitude 

Reference 

Altitude  (ft) 

10 

0.0 

6.0 

Roll  Attitude  (deg) 

10 

0.0 

0.2 

Pitch  Attitude  (deg) 

10 

0.1 

0.2 

Yaw  Attitude  (deg) 

10 

0.1 

0.2 

Engine  Thrust 

Estimated  Gross  Thrust 

10 

10.0 

300.0 

Table  A.12:  TAFA  Sensor  Model  Parameters 


A.8  Flight  Control  System  Digital  Effects 

The  flight  control  system  digital  effects  model  characterizes  the  time  delay  (i.e.  phase  lag) 
caused  by  the  computational  lag  of  the  control  laws.  The  computational  lag  is  the  time  required 
for  the  control  laws  to  generate  control  surface  commands  once  the  stick  commands  and  sensor 
feedbacks  are  available  in  the  flight  control  processor.  The  computational  lag  is  equivalent  to  a 
time  delay  which  adds  phase  lag  to  the  loop  response,  thereby  influencing  stability  and  flying 
qualities. 

The  TAFA  computational  delay  was  estimated  based  on  data  collected  under  Boeing’s 

Intelligent  Flight  Control  (IFC)  program  (NASA  funded  effort).  The  IFC  program  is  developing 

reconfigurable  flight  control  laws  which  use  system  identification  to  generate  on-line  estimates 

of  the  stability  and  control  derivatives.  This  data  is  used  to  form  linear  models  of  the  plant  at  the 

current  operating  condition.  The  plant  models  used  with  a  Riccati  equation  solver  to  generate 

linear  quadratic  control  law  gains  on-line.  The  IFC  program  is  using  a  Shark  digital  signal 

processing  (DSP)  chip  to  perform  the  system  identification,  while  the  Riccati  solution  and 

control  law  computations  are  performed  on  the  main  processor  (68040  chip).  The  Riccati 
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solution  is  performed  as  a  background  task,  which  yields  a  20  Hz  computation  rate.  The  control 
law  calculations  are  performed  at  80  Hz.  Timing  analysis  on  the  flight  hardware  has  shown 
computational  delays  for  the  calculation  of  the  control  laws  of  approximately  6  msec. 

The  Boeing  TAFA  control  laws  are  based  on  a  dynamic  inversion  control  architecture  with 
on-line  single  layer  neural  networks  to  cancel  the  error  between  the  plant  model  and  the  actual 
plant  dynamics.  System  identification  is  used  to  generate  control  derivatives  used  by  the  on  line 
control  effector  management  (ICEM)  to  generate  control  effector  commands  which  produce  the 
required  control  moments  while  optimizing  specified  performance  objectives  (load  alleviation, 

. . .).  This  approach  minimizes  the  need  for  system  identification  for  stabilizing  the  aircraft 
following  failure  or  damage.  As  a  result,  the  system  identification  is  not  time  critical  and  thus 
can  be  relegated  to  a  background  task.  Boeing’s  dynamic  inversion  control  laws  have  been 
executed  at  80  Hz  in  real  time  on  a  variety  of  flight  processors,  including  68040  and  i960  chips. 
Current  processors  (e.g.  Pentium,  Power  PC, ...)  far  exceed  the  throughput  capabilities  of  the 
68040  and  i960,  and  are  expected  to  be  able  to  handle  the  additional  overhead  of  the  on-line 
neural  networks  and  ICEM.  Timing  analysis  on  the  selected  flight  processor  will  be  required  to 
determine  whether  the  system  identification  can  be  run  as  a  background  task  on  the  main 
processor  or  whether  a  DSP  chip  will  be  required  to  host  these  algorithms.  Due  to  increased 
throughput  capability  of  current  processors,  it  is  reasonable  to  expect  that  the  reconfigurable 
flight  control  algorithms  will  exhibit  computational  delays  similar  to  those  observed  with  the  IFC 
algorithms.  A  computational  delay  of  10  msec  is  assumed  for  conservatism. 


A.9Loads  Model 

The  TAFA  loads  model  characterizes  the  shear  forces  and  toques  at  various  points  on  the 
aircraft  due  to  the  aerodynamic  loads  acting  on  the  vehicle.  The  loads  acting  on  the  vehicle  must 
remain  within  the  structural  limits  to  preserve  flight  safety  and  the  vehicle’s  fatigue  life.  The 
loads  acting  on  the  vehicle  are  a  function  of  the  flight  condition  and  control  effector  utilization. 

It  is  common  for  aircraft  control  laws  to  employ  maneuver  load  alleviation  logic,  which  utilizes 
the  control  in  a  manner  which  keeps  the  loads  within  the  structural  limits  while  achieving  the 
desired  maneuvers. 
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The  TAFA  loads  model  generates  the  shear  force,  bending  moment  and  torque  moment  at 
critical  points  on  the  aircraft  structure.  The  points  were  selected  to  be  the  canard  root,  wing  root 
and  wing  fold  as  shown  in  Figure  A.  1 1 . 


Wing  Fold 

Canard  Root 

Figure  A.11:  TAFA  Loads  Model  Definitions 

Loads  are  monitored  for  both  sides  of  the  aircraft.  The  shear  force  is  measured  positive  up  and 
expressed  in  pounds.  A  positive  bending  moment  acts  to  bend  the  tip  of  the  wing  or  canard  over 
the  top  of  the  aircraft,  while  a  positive  torque  moment  acts  to  raise  the  leading  edge  of  the  wing 
or  canard.  The  torque  and  bending  moments  are  expressed  in  inch  pounds. 

The  loads  model  was  developed  using  the  MSC/NASTRAN  doublet  lattice  aerodynamic 
method  on  a  flat  lifting  surface  model  approximating  the  TAFA  planform.  The  aerodynamic 
model  was  built  to  best  approximate  the  important  features  of  the  actual  aircraft  configuration 
given  the  constraints  of  the  doublet  lattice  aerodynamic  method.  The  model  panel  geometry  was 
designed  to  incorporate  the  canard,  leading  and  trailing  edge  flaps,  and  the  aileron.  The  layout  of 
the  wing  panel  arrangement  was  complicated  by  the  very  low  taper  ratio  and  constant  chord 
trailing  edge  control  surfaces.  The  left  half  of  the  aircraft  was  modeled,  and  symmetry  was 
assumed. 

The  aft  body  split  flaps  are  a  clam  shell  device,  with  one  panel  opening  above  the  wing  and 
the  other  below  the  wing.  The  aft  body  flaps  produce  negligible  lift  and  thus  do  not  have  a 
significant  contribution  to  the  shear  force  at  the  wing  root  or  wing  fold.  The  panels  above  and 
below  the  wing  produce  a  net  torque  and  bending  moments  near  zero  since  the  contributions  of 
the  upper  and  lower  panels  cancel  each  other.  Based  on  these  considerations,  the  aft  body  split 
flaps  were  excluded  from  the  loads  model.  The  thrust  vectoring  and  forebody  blowing  devices 
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apply  forces  directly  to  the  aircraft’s  fuselage  and  thus  do  not  contribute  to  the  aerodynamic 
loads  of  interest. 

Aerodynamic  analyses  were  conducted  to  identify  the  effect  of  the  aircraft  state  variables 
(angle  of  attack,  sideslip,  rotational  rates,  control  surface  deflections)  on  each  component  load 
for  a  variety  of  flight  conditions.  Symmetric  and  anti-symmetric  unit  trim  variable  analyses  were 
conducted  at  Mach  numbers  0.6, 0.8, 0.9, 1.1, 1.2,  and  1.4  to  produce  pressure  distributions 
across  the  planform.  These  pressures  were  then  integrated  into  canard,  wing  root,  and  wing  fold 
shear  force,  bending  moment  and  torsion  moment.  This  analysis  produced  derivatives  of  the 
desired  loads  with  respect  to  each  aircraft  state  variable.  The  loads  are  then  calculated  by 
multiplying  each  load  derivative  by  its  respective  state  variable,  and  summing  the  contributions 
across  the  state  variables. 


A.10  Atmosphere  Model 

The  TAFA  atmosphere  model  describes  the  properties  of  the  air  near  the  vehicle.  These 
properties  include  basic  atmospheric  conditions  (pressure,  temperature,  speed  of  sound,  air 
density,  and  air  coefficient  of  friction),  as  well  as  atmospheric  disturbances  (winds,  and 
turbulence).  These  models  are  described  in  the  following  sections. 


A.10.1  Atmospheric  Conditions 

The  TAFA  atmosphere  model  implements  MIL-STD-2 1 0A  and  US  Standard  Atmosphere 
1 966  Supplements  models  for  use  in  simulation  and  analysis.  The  model  allows  the  user  to 
select  from  the  various  atmosphere  models  (source  and  day  type).  In  addition,  the  user  is 
allowed  to  include  pressure  and  temperature  biases  as  a  possible  aid  in  matching  flight  test  data 
results.  The  model  produces  outside  air  temperature,  static  pressure,  air  density,  speed  of  sound, 
and  the  air  coefficient  of  friction  as  a  function  of  altitude.  The  available  atmosphere  models  are: 

MIL-STD-2 1 0A  Tropical  Day 

MIL-STD-2 1 0A  Polar  Day 

MIL-STD-2 10A  Hot  Day 
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MIL-STD-2 1 OA  Cold  Day 

U.S.  Standard  Atmosphere  Standard  Day  (45°  N  Spring/Fall) 
U.S.  Standard  Atmosphere  Hot  Day  (30°  N  July) 

U.S.  Standard  Atmosphere  Cold  Day  (75°  N  January) 


The  equations  governing  the  MIL-STD-2 10A  and  U.S.  Standard  Atmosphere  1966 
Supplements  differ,  primarily  in  the  fact  that  the  MIL-STD-2 10A  models  use  the  standard  day 
temperature  gradients  for  calculations  with  all  day  types.  Temperature  gradients  for  the  selected 
day  type  are  used  in  calculations  for  the  U.S.  Standard  Atmosphere  1966  Supplements  models. 

As  a  result,  the  U.S.  Standard  Atmosphere  1966  Supplements  models  are  implemented  via  the 
appropriate  equations,  while  the  MIL-STD-2 10A  models  are  stored  as  tables. 

The  model  produces  outside  air  temperature,  static  pressure,  air  density,  speed  of  sound,  and 
the  air  coefficient  of  friction  as  a  function  of  altitude  for  the  selected  day  type.  The  model  also 
allows  the  user  to  specify  pressure  and  temperature  biases  as  a  possible  aid  in  matching  flight  test 
data  results.  The  pressure  biases  can  be  added  in  two  ways  as  given  by  the  equation 


P  =  P_, 


(Pt  +  dPp 


+  SBias 


(A.48) 


V  ■‘o  J 

where  Pnom  is  the  nominal  pressure  profile  for  the  selected  day  type,  P0  is  the  sea  level  pressure 
for  the  selected  day  type,  dP0  is  the  sea  level  pressure  increment,  and  Pbuu  is  the  pressure  bias. 
The  sea  level  pressure  increment  produces  an  effect  which  increases  with  altitude  while 
maintaining  a  hydrostatically  baiarced  model.  The  pressure  bias  is  constant  over  all  altitudes, 
and  yields  a  model  which  is  not  hydrostatically  balanced.  The  temperature  bias  is  applied  as 
T  =  T^+Ttl„  (A.49) 


where  Tnom  is  the  nominal  temperature  profile  for  the  selected  day  type,  and  TBias  is  the 
temperature  bias.  The  temperature  bias  also  maintains  a  hydrostatically  balanced  model. 


A.10.2  Wind  and  Turbulence 

The  TAFA  simulation  is  supported  by  wind  and  turbulence  models  which  implements  the 
MIL- 1797 A  steady  state  wind,  wind  shear,  turbulence,  and  discrete  gust  characteristics  for  use 
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in  simulation  analyses.  The  model  allows  the  user  to  activate  any  combination  of  these 
atmospheric  characteristics.  The  wind  and  turbulence  models  are  described  in  detail  below. 

The  steady  state  wind  specifies  a  constant  wind  magnitude  and  direction  which  is 
independent  of  altitude.  The  wind  shear  component  models  wind  variations  as  a  function  of 
ground  clearance  altitude.  Linear,  logarithmic  and  vector  shear  models  are  available.  The  linear 
and  logarithmic  shear  models  describe  a  constant  direction  wind  whose  magnitude  varies  with 
altitude.  The  vector  shear  model  describes  a  constant  magnitude  wind  whose  direction  varies 
with  altitude.  The  turbulence  models  are  implemented  using  the  Dryden  spectra  for  linear  (axial, 
lateral,  vertical)  and  rotational  (roll,  pitch  and  yaw)  disturbances.  Light,  moderate  and  severe 
turbulence  models  are  available  for  both  high  and  low  altitude  flight.  The  discrete  gust  models 
implement  a  single  gust  profile,  triggered  at  a  specified  altitude,  specified  as  a  function  of  time, 
altitude,  or  distance  traveled. 

The  net  result  of  the  wind  and  turbulence  model  is  linear  wind  velocity  components  in 
inertial  axes,  and  rotational  gust  components  in  the  body  axes.  The  linear  wind  components  are 
combined  with  the  aircraft's  inertial  velocity  to  create  airspeed  components.  The  body  axis 
airspeed  components  are  used  to  compute  angle  of  attack  and  sideslip  angles.  The  rotational  gust 
components  are  combined  with  the  body  axis  rotational  rates  to  drive  the  damping  terms  of  the 
aerodynamics  model. 

The  implementation  of  the  angle  of  attack  and  sideslip  rate  computations  does  not  entirely 
account  for  the  turbulence  influence.  For  example,  the  angle  of  attack  rate  is  computed  as 


a  = 


uw  —  wu 
u 2  +  w2 


(A.50) 


The  wind  and  turbulence  effects  are  included  in  the  linear  velocity  terms,  but  not  the  linear 
acceleration  terms.  Analytical  computation  of  the  linear  acceleration  terms  is  not  possible  when 
turbulence  is  active  since  this  would  involve  differentiating  a  random  process.  The  influence  of 
turbulence  on  the  angle  of  attack  and  sideslip  rates  can  only  be  accurately  captured  by 
numerically  differentiating  the  corresponding  angle  of  attack  and  sideslip  angles.  An  alternate 
approach  is  to  numerically  differentiate  the  linear  velocity  terms  for  use  in  the  above  equation. 
Refined  computation  of  the  angle  of  attack  and  sideslip  rates  due  to  turbulence  is  left  for  user 
implementation  since  the  preferred  approach  may  be  application  dependent. 
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The  steady  state  wind  model  specifies  a  constant  wind  magnitude  and  direction  which  does 
not  vary  with  altitude.  The  direction  is  specified  by  an  inertial  heading  and  direction  relative  to 
the  horizontal  plane. 

The  linear  wind  shear  model  characterizes  wind  variations  as  a  linear  function  of  ground 
clearance  altitude.  The  wind  speed  varies  linearly  with  ground  clearance  altitude  with  the  wind 
maintaining  a  constant  direction.  The  wind  shear  model  contains  only  horizontal  plane  wind 
components,  with  the  vertical  wind  shear  component  being  zero.  The  linear  wind  shear  model 
parameter  definitions  are  illustrated  in  Figure  A.  12. 


Figure  A.12:  Linear  Wind  Shear  Model  Definitions 


The  loguithmic  wind  shear  model  characterizes  wind  variations  as  a  logarithmic  function  of 
ground  rl^arance  altitude.  The  wind  speed  varies  logarithmically  with  ground  clearance  altitude, 
with  the  wind  maintaining  a  constant  direction.  The  wind  shear  model  contains  only  horizontal 
plane  wind  components,  with  the  vertical  wind  shear  component  being  zero.  The  MIL_SPEC 
logarithmic  shear  model  differs  for  up-and-away  flight  and  carrier  take  off  and  landing.  The 
desired  characteristic  is  selectable  in  the  model.  The  logarithmic  wind  shear  model  has  the  form 


WS  =  WSmag 


M20/zo) 


(A.51) 


where  hcL  is  the  ground  clearance  altitude,  zo  =  2.0  for  up-and-away  and  zo  =  0. 15  for  earner 
takeoff  and  landing,  and  WS  is  the  wind  shear  magnitude.  The  up-and-away  and  carrier  wind 
shear  models  are  illustrated  in  Figure  A.  1 3 
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Figure  A.13:  Logarithmic  Wind  Shear  Model 


The  vector  wind  shear  model  describes  a  constant  magnitude  wind  whose  direction  varies 
linearly  with  ground  clearance  altitude.  The  vector  wind  shear  model  contains  only  horizontal 
plane  wind  components,  with  the  vertical  wind  shear  component  being  zero.  A  specified  change 
in  wind  direction  occurs  over  a  specified  change  in  altitude.  The  vector  wind  shear  model 
parameter  definitions  are  illustrated  in  Figure  A.  14. 


Ground  Clearance  Altitude  (ft) 


Figure  A.14:  Vector  Wind  Shear  Model 
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The  continuous  turbulence  models  are  implemented  using  the  Dryden  spectra.  Light, 
moderate  and  severe  turbulence  models  are  available  for  both  high  and  low  altitude  flight  The 
spectra  are  implemented  using  weighting  filters  driven  by  Gaussian  white  noise  sources  to 
produce  the  required  correlation  characteristics.  Linear  (axial,  lateral,  and  vertical)  and 
rotational  (roll,  pitch,  and  yaw)  turbulence  models  are  implemented.  The  linear  turbulence  is 
directed  along  the  horizontal  flight  path  angle,  while  the  rotational  turbulence  spectra  are 
implemented  in  body  axes.  Gust  multipliers  are  provided  to  scale  (or  turn  off)  the  turbulence  in 
selected  axes. 

The  Dryden  spectra  and  continuous  domain  weighting  filter  implementations  are  summarized 
in  Figure  A.  15. 
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Figure  A.15:  Dryden  Turbulence  Spectra  and  Continuous  Weighting  Filter 
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The  turbulence  model  parameters  for  light,  moderate,  and  severe  turbulence  are  provided  in 
Table  A.  13. 
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Condition 

RMS  Intensities  ~  ft/sec 

(au,  Gv>  <*w) 

Scale  Lengths  ~  ft 
(Lu,  Lv,  Lw) 

Altitude  >  1000  ft 

>  ftinction  of 

altitude  and  turbulence 

severity  shown  in  Figure 

A.16 

Lu  = 1750 

Lv  =  Lw  =  875 

Altitude  <  1000  ft 

(Carrier 

Environment) 

Gu=1.77,  Gv=1-69,  G\y=1.06 

Lu=Lv=Lw=100 

Altitude  <  1000  ft 

(Land 

Environment) 

aw=0.1*U20 

au=<Jw/(0. 1 77+0.000823h)°- 

4 

Lu=h/(0. 1 77+0.000823h)  1  1 
Lv=Lu/2 

Lw=h/2 

U20  =  15  knots  (light  turb.),  30  knots  (moderate  turb.),  45  knots  (severe  turb.) 


Table  A.13:  Turbulence  Model  Parameters 
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Figure  A.16:  Medium  /  High  Altitude  Turbulence  Intensities 
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The  discrete  gust  models  implement  a  single  gust  profile,  triggered  at  a  specified  altitude,  and 
specified  as  a  function  of  time,  altitude,  or  distance  traveled.  The  time  dependent  gusts  can  be 
used  to  implement  user  specified  gust  profiles  in  the  axial,  lateral,  vertical  and  roll  axes.  The 
altitude  and  distance  dependent  gusts  exhibit  a  one  minus  cosine  characteristic  and  can  be 
employed  in  one  of  the  four  axes  (axial,  lateral,  vertical  and  roll).  The  user  is  allowed  to  specify 
the  frequency  of  the  gust  profile,  the  maximum  gust  magnitude,  and  distance  between  gusts. 

The  discrete  gust  as  a  function  of  time  model  implements  a  single  gust  profile,  triggered  at  a 
specified  altitude.  If  a  gust  profile  is  not  provided,  the  MIL-SPEC  discrete  gust  profile  of  Figure 
A.  17  is  applied  in  the  lateral  and  vertical  directions.  The  gust  profile  of  Figure  A.  17  represents 
an  actual  wind  gust  encountered  during  air  vehicle  flight  testing. 


o 


Time  From  Start  of  Gust  -  Sec 
Figure  A.17:  MIL-SPEC  Discrete  Gust  Profile 


The  discrete  gust  as  a  function  of  altitude  or  distance  models  implement  a  single  gust  profile, 
triggered  at  a  specified  altitude.  These  gusts  exhibit  a  one  minus  cosine  characteristic  and  can  be 
employed  in  one  of  the  four  axes  (axial,  lateral,  vertical  and  roll).  The  user  is  allowed  to  specify 
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the  frequency  of  the  gust  profile,  the  maximum  gust  magnitude,  and  distance  between  gusts. 
Definitions  of  these  parameters  are  shown  in  Figure  A.  1 8. 


Figure  A.18:  One  Minus  Cosine  Discrete  Gust  Model 


